Index: /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog_nogcm_pluto.F =================================================================== --- /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog_nogcm_pluto.F (revision 3229) +++ /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog_nogcm_pluto.F (revision 3229) @@ -0,0 +1,1048 @@ +! +! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $ +! +c +c + SUBROUTINE leapfrog_nogcm_pluto(ucov,vcov,teta,ps,masse,phis,q, + & time_0) + + +cIM : pour sortir les param. du modele dans un fis. netcdf 110106 +#ifdef CPP_IOIPSL + use IOIPSL +#endif + USE infotrac, ONLY: nqtot,ok_iso_verif,tname + USE guide_mod, ONLY : guide_main + USE write_field, ONLY: writefield + USE control_mod, ONLY: planet_type,nday,day_step,iperiod,iphysiq, + & less1day,fractday,ndynstep,iconser, + & dissip_period,offline,ip_ebil_dyn, + & ok_dynzon,periodav,ok_dyn_ave,iecri, + & ok_dyn_ins,output_grads_dyn + use exner_hyb_m, only: exner_hyb + use exner_milieu_m, only: exner_milieu + use cpdet_mod, only: cpdet,tpot2t,t2tpot + use sponge_mod, only: callsponge,mode_sponge,sponge + use comuforc_h + USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs, + & aps,bps,presnivs,pseudoalt,preff,scaleheight + USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss, + . cpp,ihf,iflag_top_bound,pi,kappa,r + USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys, + . statcl,conser,purmats,tidal,ok_strato + USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref, + . start_time,dt + + + + + 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 pour l'advection de +c q , en modifiant iadv dans traceur.def (10/02) . +c +c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99) +c Pour Van-Leer iadv=10 +c +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +#include "dimensions.h" +#include "paramet.h" +#include "comdissnew.h" +#include "comgeom.h" +!#include "com_io_dyn.h" +#include "iniprint.h" +#include "academic.h" + +! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique +! #include "clesphys.h" + + REAL,INTENT(IN) :: time_0 ! not used + +c dynamical variables: + REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind + REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind + REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature + REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa) + REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass + REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface + REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers + + REAL p (ip1jmp1,llmp1 ) ! interlayer pressure + REAL pks(ip1jmp1) ! exner at the surface + REAL pk(ip1jmp1,llm) ! exner at mid-layer + REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer + REAL phi(ip1jmp1,llm) ! geopotential + REAL w(ip1jmp1,llm) ! vertical velocity + REAL kpd(ip1jmp1) ! TB15 = exp (-z/H) + +! ADAPTATION GCM POUR CP(T) + REAL temp(ip1jmp1,llm) ! temperature + REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk + +! real zqmin,zqmax + + +! ED18 nogcm + REAL tau_ps + REAL tau_n2 + REAL tau_teta + REAL tetadpmean + REAL tetadp(ip1jmp1,llm) + REAL dpmean + REAL dp(ip1jmp1,llm) + REAL tetamean + real globaverage2d +! LOGICAL firstcall_globaverage2d + + +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 en */s + REAL dv(ip1jm,llm),du(ip1jmp1,llm) + REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot) + +c tendances de la dissipation en */s + REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm) + REAL dtetadis(ip1jmp1,llm) + +c tendances de la couche superieure */s +c REAL dvtop(ip1jm,llm) + REAL dutop(ip1jmp1,llm) +c REAL dtetatop(ip1jmp1,llm) +c REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1) + +c TITAN : tendances due au forces de marees */s + REAL dvtidal(ip1jm,llm),dutidal(ip1jmp1,llm) + +c tendances physiques */s + REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm) + REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1) + +c variables pour le fichier histoire +! REAL dtav ! intervalle de temps elementaire + + REAL tppn(iim),tpps(iim),tpn,tps +c + INTEGER itau,itaufinp1,iav +! INTEGER iday ! jour julien + REAL time + + REAL SSUM +! REAL finvmaold(ip1jmp1,llm) + +cym LOGICAL lafin + LOGICAL :: lafin=.false. + INTEGER ij,iq,l +! INTEGER ik + +! real time_step, t_wrt, t_ops + + REAL rdaym_ini +! jD_cur: jour julien courant +! jH_cur: heure julienne courante + REAL :: jD_cur, jH_cur +! INTEGER :: an, mois, jour +! REAL :: secondes + + LOGICAL first,callinigrads +cIM : pour sortir les param. du modele dans un fis. netcdf 110106 + save first + data first/.true./ +! real dt_cum +! character*10 infile +! integer zan, tau0, thoriid +! integer nid_ctesGCM +! save nid_ctesGCM +! real degres +! real rlong(iip1), rlatg(jjp1) +! real zx_tmp_2d(iip1,jjp1) +! integer ndex2d(iip1*jjp1) + logical ok_sync + parameter (ok_sync = .true.) + logical physics + + data callinigrads/.true./ + character*10 string10 + +! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) + REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale + + REAL psmean ! pression moyenne + REAL pqn2mean ! moyenne globale ps*qco + REAL p0 ! pression de reference + REAL p00d ! globalaverage(kpd) + REAL qmean_n2,qmean_n2_vert ! mass mean mixing ratio vap n2 + REAL pqn2(ip1jmp1) ! average n2 mass index : ps*q_n2 + REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4 + +c+jld variables test conservation energie + REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm) +C Tendance de la temp. potentiel d (theta)/ d t due a la +C tansformation d'energie cinetique en energie thermique +C cree par la dissipation + REAL dtetaecdt(ip1jmp1,llm) + REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm) + REAL vnat(ip1jm,llm),unat(ip1jmp1,llm) +! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec + CHARACTER*15 ztit +!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag. +!IM SAVE ip_ebil_dyn +!IM DATA ip_ebil_dyn/0/ +c-jld + +! integer :: itau_w ! for write_paramLMDZ_dyn.h + +! character*80 dynhist_file, dynhistave_file + character(len=*),parameter :: modname="leapfrog" + character*80 abort_message + + logical dissip_conservative + save dissip_conservative + data dissip_conservative/.true./ + + INTEGER testita + PARAMETER (testita = 9) + + logical , parameter :: flag_verif = .false. + + ! for CP(T) + real :: dtec + real :: ztetaec(ip1jmp1,llm) + + ! TEMP : diagnostic mass + real :: n2mass(iip1,jjp1) + real :: n2ice_ij(iip1,jjp1) + integer,save :: igcm_n2=0 ! index of CO2 tracer (if any) + integer :: i,j,ig + integer, parameter :: ngrid = 2+(jjm-1)*iim + ! local mass + real :: mq(ip1jmp1,llm) + + if (nday>=0) then + itaufin = nday*day_step + else + ! to run a given (-nday) number of dynamical steps + itaufin = -nday + endif + if (less1day) then +c MODIF VENUS: to run less than one day: + itaufin = int(fractday*day_step) + endif + if (ndynstep.gt.0) then + ! running a given number of dynamical steps + itaufin=ndynstep + endif + itaufinp1 = itaufin +1 + + + + +c INITIALISATIONS + dudis(:,:) =0. + dvdis(:,:) =0. + dtetadis(:,:)=0. + dutop(:,:) =0. +c dvtop(:,:) =0. +c dtetatop(:,:)=0. +c dqtop(:,:,:) =0. +c dptop(:) =0. + dufi(:,:) =0. + dvfi(:,:) =0. + dtetafi(:,:)=0. + dqfi(:,:,:) =0. + dpfi(:) =0. + + itau = 0 + physics=.true. + if (iflag_phys==0.or.iflag_phys==2) physics=.false. + +! ED18 nogcm +! firstcall_globaverage2d = .true. + +c iday = day_ini+itau/day_step +c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0 +c IF(time.GT.1.) THEN +c time = time-1. +c iday = iday+1 +c ENDIF + + +c----------------------------------------------------------------------- +c On initialise la pression et la fonction d'Exner : +c -------------------------------------------------- + + dq(:,:,:)=0. + CALL pression ( ip1jmp1, ap, bp, ps, p ) + if (pressure_exner) then + CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf ) + else + CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf ) + endif + +c------------------ +c TEST PK MONOTONE +c------------------ + write(*,*) "Test PK" + do ij=1,ip1jmp1 + do l=2,llm +! PRINT*,'pk(ij,l) = ',pk(ij,l) +! PRINT*,'pk(ij,l-1) = ',pk(ij,l-1) + if(pk(ij,l).gt.pk(ij,l-1)) then +c write(*,*) ij,l,pk(ij,l) + abort_message = 'PK non strictement decroissante' + call abort_gcm(modname,abort_message,1) +c write(*,*) "ATTENTION, Test PK deconnecte..." + endif + enddo + enddo + write(*,*) "Fin Test PK" +c stop +c------------------ +c Preparing mixing of pressure and tracers in nogcm +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*200.)) + ENDDO + p00d=globaverage2d(kpd) ! mean pres at ref level + tau_ps = 1. ! constante de rappel for pressure (s) + tau_n2 = 1.E5 !E5 ! constante de rappel for mix ratio qn2 (s) + tau_teta = 1.E7 !constante de rappel for potentiel temperature + + PRINT*,'igcm_n2 = ',igcm_n2 + do iq=1,nqtot + if (tname(iq)=="n2") then + igcm_n2=iq + exit + endif + enddo + +c----------------------------------------------------------------------- +c Debut de l'integration temporelle: +c ---------------------------------- + + 1 CONTINUE ! Matsuno Forward step begins here + +c date: (NB: date remains unchanged for Backward step) +c ----- + + jD_cur = jD_ref + day_ini - day_ref + & + & (itau+1)/day_step + jH_cur = jH_ref + start_time + & + & mod(itau+1,day_step)/float(day_step) + jD_cur = jD_cur + int(jH_cur) + jH_cur = jH_cur - int(jH_cur) + +c + +! Save fields obtained at previous time step as '...m1' + 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 + + 2 CONTINUE ! Matsuno backward or leapfrog step begins here + +c----------------------------------------------------------------------- + +c date: (NB: only leapfrog step requires recomputing date) +c ----- + + IF (leapf) THEN + jD_cur = jD_ref + day_ini - day_ref + + & (itau+1)/day_step + jH_cur = jH_ref + start_time + + & mod(itau+1,day_step)/float(day_step) + jD_cur = jD_cur + int(jH_cur) + jH_cur = jH_cur - int(jH_cur) + ENDIF + + +c gestion des appels de la physique et des dissipations: +c ------------------------------------------------------ +c +c ... P.Le Van ( 6/02/95 ) .... + +! ED18: suppression des mentions de la variable apdiss dans le cas +! 'nogcm' + + apphys = .FALSE. + statcl = .FALSE. + conser = .FALSE. + + + IF( purmats ) THEN + ! Purely Matsuno time stepping + IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE. + + + IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward + s .and. physics ) apphys = .TRUE. + ELSE + ! Leapfrog/Matsuno time stepping + IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE. + + + IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE. + END IF + +! Ehouarn: for Shallow Water case (ie: 1 vertical layer), +! supress dissipation step + + + +c----------------------------------------------------------------------- +c calcul des tendances dynamiques: +c -------------------------------- + +! ED18: suppression de l'onglet pour le cas nogcm + + + dv(:,:) = 0. + du(:,:) = 0. + dteta(:,:) = 0. + dq(:,:,:) = 0. + + + +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 ) + if (pressure_exner) then + CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf ) + else + CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf ) + endif + +! Compute geopotential (physics might need it) + CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) + + jD_cur = jD_ref + day_ini - day_ref + & + & (itau+1)/day_step + + ! AS: we make jD_cur to be pday + jD_cur = int(day_ini + itau/day_step) + +! print*,'itau =',itau +! print*,'day_step =',day_step +! print*,'itau/day_step =',itau/day_step + + + jH_cur = jH_ref + start_time + & + & mod(itau+1,day_step)/float(day_step) + jH_cur = jH_ref + start_time + & + & mod(itau,day_step)/float(day_step) + jD_cur = jD_cur + int(jH_cur) + jH_cur = jH_cur - int(jH_cur) + + + +! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur +! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes) +! write(lunout,*)'current date = ',an, mois, jour, secondes + +c rajout debug +c lafin = .true. + + +c Interface avec les routines de phylmd (phymars ... ) +c ----------------------------------------------------- + +c+jld + +c Diagnostique de conservation de l'Energie : initialisation + IF (ip_ebil_dyn.ge.1 ) THEN + ztit='bil dyn' +! Ehouarn: be careful, diagedyn is Earth-specific! + IF (planet_type.eq."earth") THEN + CALL diagedyn(ztit,2,1,1,dtphys + & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2)) + ENDIF + ENDIF ! of IF (ip_ebil_dyn.ge.1 ) +c-jld +#ifdef CPP_IOIPSL +cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ +cIM uncomment next 6 lines to get some parameters for LMDZ dynamics +c IF (first) THEN +c first=.false. +c#include "ini_paramLMDZ_dyn.h" +c ENDIF +c +c#include "write_paramLMDZ_dyn.h" +c +#endif +! #endif of #ifdef CPP_IOIPSL + +c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/))) +! print*,'---------LEAPFROG---------------' +! print*,'' +! print*,'> AVANT CALFIS' +! print*,'' +! print*,'teta(3049,:) = ',teta(3049,:) +! print*,'' +! print*,'dteta(3049,:) = ',dteta(3049,:) +! print*,'' +! print*,'dtetafi(3049,:) = ',dtetafi(3049,:) +! print*,'' + + + CALL calfis( lafin , jD_cur, jH_cur, + $ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi , + $ du,dv,dteta,dq, + $ flxw, + $ dufi,dvfi,dtetafi,dqfi,dpfi ) + +c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/))) +c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/))) +c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/))) +! print*,'> APRES CALFIS (AVANT ADDFI)' +! print*,'' +! print*,'teta(3049,:) = ',teta(3049,:) +! print*,'' +! print*,'dteta(3049,:) = ',dteta(3049,:) +! print*,'' +! print*,'dtetafi(3049,:) = ',dtetafi(3049,:) +! print*,'' + + +c ajout des tendances physiques +c ------------------------------ + + CALL addfi( dtphys, leapf, forward , + $ ucov, vcov, teta , q ,ps , + $ dufi, dvfi, dtetafi , dqfi ,dpfi ) + +! print*,'> APRES ADDFI' +! print*,'' +! print*,'teta(3049,:) = ',teta(3049,:) +! print*,'' +! print*,'dteta(3049,:) = ',dteta(3049,:) +! print*,'' +! print*,'dtetafi(3049,:) = ',dtetafi(3049,:) +! print*,'' + + CALL pression (ip1jmp1,ap,bp,ps,p) + CALL massdair(p,masse) + + ! Mass of tracers in each mess (kg) before Horizontal mixing of pressure +c ------------------------------------------------------------------- + + DO l=1, llm + DO ij=1,ip1jmp1 + mq(ij,l) = masse(ij,l)*q(ij,l,igcm_n2) + ENDDO + ENDDO + +c Horizontal mixing of pressure +c ------------------------------ + ! Rappel newtonien vers psmean + psmean= globaverage2d(ps) ! mean pressure + p0=psmean/p00d + DO ij=1,ip1jmp1 + oldps(ij)=ps(ij) + ps(ij)=ps(ij) +(p0*kpd(ij) + & -ps(ij))*(1-exp(-dtphys/tau_ps)) + ENDDO + + write(72,*) 'psmean = ',psmean ! mean pressure redistributed + + ! security check: 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*kpd(ij) + !write(*,*) 'tau_ps = ',tau_ps + !STOP + ps(ij)=0.0000001*kpd(ij)/p00d + ENDIF + ENDDO +!*********************** +! Correction on teta due to surface pressure changes + DO l = 1,llm + DO ij = 1,ip1jmp1 + teta(ij,l)= teta(ij,l)*(ps(ij)/oldps(ij))**kappa + ENDDO + ENDDO +!*********************** + + + +! ! update pressure and update p and pk +! DO ij=1,ip1jmp1 +! ps(ij) = ps(ij) + dpfi(ij)*dtphys +! ENDDO + CALL pression (ip1jmp1,ap,bp,ps,p) + CALL massdair(p,masse) + if (pressure_exner) then + CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf ) + else + CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf ) + endif + + ! Update tracers after Horizontal mixing of pressure to ! conserve tracer mass +c ------------------------------------------------------------------- + DO l=1, llm + DO ij=1,ip1jmp1 + q(ij,l,igcm_n2) = mq(ij,l)/ masse(ij,l) + ENDDO + ENDDO + + + +c Horizontal mixing of pressure +c ------------------------------ + ! Rappel newtonien vers psmean + psmean= globaverage2d(ps) ! mean pressure +! ! increment q_n2 with physical tendancy +! IF (igcm_n2.ne.0) then +! DO l=1, llm +! DO ij=1,ip1jmp1 +! q(ij,l,igcm_n2)=q(ij,l,igcm_n2)+ +! & dqfi(ij,l,igcm_n2)*dtphys +! ENDDO +! ENDDO +! ENDIF + +c Mixing CO2 vertically +c -------------------------- + if (igcm_n2.ne.0) then + DO ij=1,ip1jmp1 + qmean_n2_vert=0. + DO l=1, llm + qmean_n2_vert= qmean_n2_vert + & + q(ij,l,igcm_n2)*( p(ij,l) - p(ij,l+1)) + END DO + qmean_n2_vert= qmean_n2_vert/ps(ij) + DO l=1, llm + q(ij,l,igcm_n2)= qmean_n2_vert + END DO + END DO + end if + + + +c Horizontal mixing of pressure +c ------------------------------ + ! Rappel newtonien vers psmean + psmean= globaverage2d(ps) ! mean pressure + +c Horizontal mixing tracers, and temperature (replace dynamics in nogcm) +c --------------------------------------------------------------------------------- + + ! Simulate redistribution by dynamics for qn2 + if (igcm_n2.ne.0) then + + DO ij=1,ip1jmp1 + pqn2(ij)= ps(ij) * q(ij,1,igcm_n2) + ENDDO + pqn2mean=globaverage2d(pqn2) + + ! Rappel newtonien vers qn2_mean + qmean_n2= pqn2mean / psmean + + DO ij=1,ip1jmp1 + q(ij,1,igcm_n2)=q(ij,1,igcm_n2)+ + & (qmean_n2-q(ij,1,igcm_n2))* + & (1.-exp(-dtphys/tau_n2)) + ENDDO + + DO l=2, llm + DO ij=1,ip1jmp1 + q(ij,l,igcm_n2)=q(ij,1,igcm_n2) + END DO + END DO + +! TEMPORAIRE (ED) +! PRINT*,'psmean = ',psmean +! PRINT*,'qmean_n2 = ',qmean_n2 +! PRINT*,'pqn2mean = ',pqn2mean +! PRINT*,'q(50,1,igcm_n2) = ',q(50,1,igcm_n2) +! PRINT*,'q(50,2,igcm_n2) = ',q(50,2,igcm_n2) +! PRINT*,'q(50,3,igcm_n2) = ',q(50,3,igcm_n2) + + endif ! igcm_n2.ne.0 + + +! ******************************************************** + +c Horizontal mixing of pressure +c ------------------------------ + ! Rappel newtonien vers psmean + psmean= globaverage2d(ps) ! mean pressure + + +c Mixing Temperature horizontally +c ------------------------------- +! initialize variables that will be averaged + DO l=1,llm + DO ij=1,ip1jmp1 + dp(ij,l) = p(ij,l) - p(ij,l+1) + tetadp(ij,l) = teta(ij,l)*dp(ij,l) + ENDDO + ENDDO + + DO l=1,llm + tetadpmean = globaverage2d(tetadp(:,l)) + dpmean = globaverage2d(dp(:,l)) + tetamean = tetadpmean / dpmean + DO ij=1,ip1jmp1 + teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) * + & (1 - exp(-dtphys/tau_teta)) + ENDDO + ENDDO + + + ENDIF ! of IF( apphys ) + + +c ******************************************************************** +c ******************************************************************** +c .... fin de l'integration dynamique et physique pour le pas itau .. +c ******************************************************************** +c ******************************************************************** + + + IF ( .NOT.purmats ) THEN +c ........................................................ +c .............. schema matsuno + leapfrog .............. +c ........................................................ + + IF(forward. OR. leapf) THEN + itau= itau + 1 +c iday= day_ini+itau/day_step +c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0 +c IF(time.GT.1.) THEN +c time = time-1. +c iday = iday+1 +c ENDIF + ENDIF + + + IF( itau. EQ. itaufinp1 ) then + if (flag_verif) then + write(79,*) 'ucov',ucov + write(80,*) 'vcov',vcov + write(81,*) 'teta',teta + write(82,*) 'ps',ps + write(83,*) 'q',q + WRITE(85,*) 'q1 = ',q(:,:,1) + WRITE(86,*) 'q3 = ',q(:,:,3) + endif + + 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 + +! ! Ehouarn: re-compute geopotential for outputs + CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) + + IF (ok_dynzon) THEN +#ifdef CPP_IOIPSL +c les traceurs ne sont pas sortis, trop lourd. +c Peut changer eventuellement si besoin. + CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav, + & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov, + & du,dudis,dutop,dufi) +#endif + END IF + IF (ok_dyn_ave) THEN +#ifdef CPP_IOIPSL + CALL writedynav(itau,vcov, + & ucov,teta,pk,phi,q,masse,ps,phis) +#endif + ENDIF + + ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin)) + + if (ok_iso_verif) then + call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584') + endif !if (ok_iso_verif) then + +c----------------------------------------------------------------------- +c ecriture de la bande histoire: +c ------------------------------ + + IF( MOD(itau,iecri).EQ.0) THEN + ! Ehouarn: output only during LF or Backward Matsuno + if (leapf.or.(.not.leapf.and.(.not.forward))) then +! ADAPTATION GCM POUR CP(T) + call tpot2t(ijp1llm,teta,temp,pk) + tsurpk = cpp*temp/pk + CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi) + unat=0. + do l=1,llm + unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm) + vnat(:,l)=vcov(:,l)/cv(:) + enddo +#ifdef CPP_IOIPSL + if (ok_dyn_ins) then +! write(lunout,*) "leapfrog: call writehist, itau=",itau + CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis) +! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/))) +! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/))) +! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/))) +! call WriteField('ps',reshape(ps,(/iip1,jmp1/))) +! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/))) + endif ! of if (ok_dyn_ins) +#endif +! For some Grads outputs of fields + if (output_grads_dyn) then +#include "write_grads_dyn.h" + endif + endif ! of if (leapf.or.(.not.leapf.and.(.not.forward))) + ENDIF ! of IF(MOD(itau,iecri).EQ.0) + + IF(itau.EQ.itaufin) THEN + + CALL dynredem1("restart.nc",start_time, + & vcov,ucov,teta,q,masse,ps) + CLOSE(99) + !!! Ehouarn: Why not stop here and now? + ENDIF ! of IF (itau.EQ.itaufin) + +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 ! of IF (forward) + ELSE + +c ...... pas leapfrog ..... + + leapf = .TRUE. + dt = 2.*dtvr + GO TO 2 + END IF ! of IF (MOD(itau,iperiod).EQ.0) + ! ELSEIF (MOD(itau-1,iperiod).EQ.0) + + ELSE ! of IF (.not.purmats) + + if (ok_iso_verif) then + call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664') + endif !if (ok_iso_verif) then + +c ........................................................ +c .............. schema matsuno ............... +c ........................................................ + IF( forward ) THEN + + itau = itau + 1 +c iday = day_ini+itau/day_step +c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0 +c +c IF(time.GT.1.) THEN +c time = time-1. +c iday = iday+1 +c ENDIF + + forward = .FALSE. + IF( itau. EQ. itaufinp1 ) then + abort_message = 'Simulation finished' + call abort_gcm(modname,abort_message,0) + ENDIF + GO TO 2 + + ELSE ! of IF(forward) i.e. backward step + + if (ok_iso_verif) then + call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698') + endif !if (ok_iso_verif) then + + IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN + IF(itau.EQ.itaufin) THEN + iav=1 + ELSE + iav=0 + ENDIF + +! ! Ehouarn: re-compute geopotential for outputs + CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) + + IF (ok_dynzon) THEN +#ifdef CPP_IOIPSL +c les traceurs ne sont pas sortis, trop lourd. +c Peut changer eventuellement si besoin. + CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav, + & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov, + & du,dudis,dutop,dufi) +#endif + ENDIF + IF (ok_dyn_ave) THEN +#ifdef CPP_IOIPSL + CALL writedynav(itau,vcov, + & ucov,teta,pk,phi,q,masse,ps,phis) +#endif + ENDIF + + ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) + + IF(MOD(itau,iecri ).EQ.0) THEN +c IF(MOD(itau,iecri*day_step).EQ.0) THEN +! ADAPTATION GCM POUR CP(T) + call tpot2t(ijp1llm,teta,temp,pk) + tsurpk = cpp*temp/pk + CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi) + unat=0. + do l=1,llm + unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm) + vnat(:,l)=vcov(:,l)/cv(:) + enddo +#ifdef CPP_IOIPSL + if (ok_dyn_ins) then +! write(lunout,*) "leapfrog: call writehist (b)", +! & itau,iecri + CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis) + endif ! of if (ok_dyn_ins) +#endif +! For some Grads outputs + if (output_grads_dyn) then +#include "write_grads_dyn.h" + endif + + ENDIF ! of IF(MOD(itau,iecri ).EQ.0) + + IF(itau.EQ.itaufin) THEN + CALL dynredem1("restart.nc",start_time, + & vcov,ucov,teta,q,masse,ps) + ENDIF ! of IF(itau.EQ.itaufin) + + forward = .TRUE. + GO TO 1 + + ENDIF ! of IF (forward) + + END IF ! of IF(.not.purmats) + + STOP + END + +! ************************************************************ +! globalaverage VERSION PLuto +! ************************************************************ + 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=2,jjp1-1 ! lat + DO i = 1, iim ! lon + airetot = airetot + aire(i,j) + END DO + END DO + DO i=1,iim + airetot=airetot + aire(i,1)+aire(i,jjp1) + ENDDO + 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 + + DO i=1,iim + globaverage2d = globaverage2d + var(i,1)*aire(i,1) + globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1) + ENDDO + + globaverage2d = globaverage2d/airetot + return + end + Index: /trunk/LMDZ.COMMON/libf/dyn3d/nogcm_pluto.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/dyn3d/nogcm_pluto.F90 (revision 3229) +++ /trunk/LMDZ.COMMON/libf/dyn3d/nogcm_pluto.F90 (revision 3229) @@ -0,0 +1,427 @@ +! +! $Id: gcm.F 1446 2010-10-22 09:27:25Z emillour $ +! +! En cours de modification par T. Bertrand +! +PROGRAM nogcm_pluto + +#ifdef CPP_IOIPSL + USE IOIPSL +#else + ! if not using IOIPSL, we still need to use (a local version of) getin + USE ioipsl_getincom +#endif + + +#ifdef CPP_XIOS + ! ug Pour les sorties XIOS + USE wxios +#endif + + USE filtreg_mod + USE infotrac + USE control_mod, only: planet_type,nday,day_step,iperiod,iphysiq, & + raz_date,anneeref,starttime,dayref, & + ok_dyn_ins,ok_dyn_ave,iecri,periodav, & + less1day,fractday,ndynstep,nsplit_phys + USE mod_const_mpi, ONLY: COMM_LMDZ + use cpdet_mod, only: ini_cpdet + USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref, & + itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique +! A nettoyer. On ne veut qu'une ou deux routines d'interface +! dynamique -> physique pour l'initialisation +! Ehouarn: the following are needed with (parallel) physics: +#ifdef CPP_PHYS + USE iniphysiq_mod, ONLY: iniphysiq +#endif + + USE comconst_mod, ONLY: daysec,dtvr,dtphys,rad,g,r,cpp + USE logic_mod, ONLY: read_start,iflag_phys,ok_guide,ecripar + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + IMPLICIT NONE + + ! ...... Version du 10/01/98 .......... + + ! avec coordonnees verticales hybrides + ! avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 ) + + !======================================================================= + ! + ! Auteur: P. Le Van /L. Fairhead/F.Hourdin + ! ------- + ! + ! Objet: + ! ------ + ! + ! GCM LMD nouvelle grille + ! + !======================================================================= + ! + ! ... Dans inigeom , nouveaux calculs pour les elongations cu , cv + ! et possibilite d'appeler une fonction f(y) a derivee tangente + ! hyperbolique a la place de la fonction a derivee sinusoidale. + ! ... Possibilite de choisir le schema pour l'advection de + ! q , en modifiant iadv dans traceur.def (MAF,10/02) . + ! + ! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99) + ! Pour Van-Leer iadv=10 + ! + !----------------------------------------------------------------------- + ! Declarations: + ! ------------- + + include "dimensions.h" + include "paramet.h" + include "comdissnew.h" + include "comgeom.h" +!!!!!!!!!!!#include "control.h" +!#include "com_io_dyn.h" + include "iniprint.h" + include "tracstoke.h" + + REAL zdtvr + + ! variables dynamiques + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants + REAL teta(ip1jmp1,llm) ! temperature potentielle + REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes + REAL ps(ip1jmp1) ! pression au sol + REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.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 + + ! variables dynamiques intermediaire pour le transport + + ! variables pour le fichier histoire + REAL dtav ! intervalle de temps elementaire + + REAL time_0 + + LOGICAL lafin + INTEGER ij,iq,l,i,j + + + real time_step, t_wrt, t_ops + + LOGICAL first + + ! LOGICAL call_iniphys + ! data call_iniphys/.true./ + + !+jld variables test conservation energie + ! REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm) + ! Tendance de la temp. potentiel d (theta)/ d t due a la + ! tansformation d'energie cinetique en energie thermique + ! cree par la dissipation + REAL dhecdt(ip1jmp1,llm) + ! REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm) + ! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec + CHARACTER (len=15) :: ztit + !-jld + + + character (len=80) :: dynhist_file, dynhistave_file + character (len=20) :: modname + character (len=80) :: abort_message + ! locales pour gestion du temps + INTEGER :: an, mois, jour + REAL :: heure + logical use_filtre_fft + +!----------------------------------------------------------------------- +! Initialisations: +! ---------------- + + abort_message = 'last timestep reached' + modname = 'gcm' + lafin = .FALSE. + dynhist_file = 'dyn_hist.nc' + dynhistave_file = 'dyn_hist_ave.nc' + + + +!---------------------------------------------------------------------- +! lecture des fichiers gcm.def ou run.def +! --------------------------------------- +! + CALL conf_gcm( 99, .TRUE. ) + if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", & + "iphysiq must be a multiple of iperiod", 1) + + use_filtre_fft=.FALSE. + CALL getin('use_filtre_fft',use_filtre_fft) + IF (use_filtre_fft) call abort_gcm("gcm",'FFT filter is not available in the ' & + // 'sequential version of the dynamics.', 1) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Initialisation de XIOS +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#ifdef CPP_XIOS + CALL wxios_init("LMDZ") +#endif + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! FH 2008/05/02 +! A nettoyer. On ne veut qu'une ou deux routines d'interface +! dynamique -> physique pour l'initialisation +!#ifdef CPP_PHYS +! CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) +!! call initcomgeomphy ! now done in iniphysiq +!#endif +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! +! Initialisations pour Cp(T) Venus + call ini_cpdet +! +!----------------------------------------------------------------------- +! Choix du calendrier +! ------------------- + +! calend = 'earth_365d' + +#ifdef CPP_IOIPSL + if (calend == 'earth_360d') then + call ioconf_calendar('360d') + write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an' + else if (calend == 'earth_365d') then + call ioconf_calendar('noleap') + write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an' + else if (calend == 'earth_366d') then + call ioconf_calendar('gregorian') + write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile' + else + abort_message = 'Mauvais choix de calendrier' + call abort_gcm(modname,abort_message,1) + endif +#endif +!----------------------------------------------------------------------- + ! + ! + !------------------------------------ + ! Initialisation partie parallele + !------------------------------------ + + ! + ! + !----------------------------------------------------------------------- + ! Initialisation des traceurs + ! --------------------------- + ! Choix du nombre de traceurs et du schema pour l'advection + ! dans fichier traceur.def, par default ou via INCA + call infotrac_init + + ! Allocation de la tableau q : champs advectes + allocate(q(ip1jmp1,llm,nqtot)) + + !----------------------------------------------------------------------- + ! Lecture de l'etat initial : + ! --------------------------- + + ! lecture du fichier start.nc + if (read_start) then + ! we still need to run iniacademic to initialize some + ! constants & fields, if we run the 'newtonian' or 'SW' cases: + if (iflag_phys.ne.1) then + CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) + endif + + CALL dynetat0("start.nc",vcov,ucov, & + teta,q,masse,ps,phis, time_0) + + ! write(73,*) 'ucov',ucov + ! write(74,*) 'vcov',vcov + ! write(75,*) 'teta',teta + ! write(76,*) 'ps',ps + ! write(77,*) 'q',q + + endif ! of if (read_start) + + + ! le cas echeant, creation d un etat initial + IF (prt_level > 9) WRITE(lunout,*) & + 'NOGCM: AVANT iniacademic AVANT AVANT AVANT AVANT' + if (.not.read_start) then + CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) + endif + + + !----------------------------------------------------------------------- + ! Lecture des parametres de controle pour la simulation : + ! ------------------------------------------------------- + ! on recalcule eventuellement le pas de temps + + IF(MOD(day_step,iperiod).NE.0) THEN + abort_message = & + 'Il faut choisir un nb de pas par jour multiple de iperiod' + call abort_gcm(modname,abort_message,1) + ENDIF + +! IF(MOD(day_step,iphysiq).NE.0) THEN +! abort_message = & +! 'Il faut choisir un nb de pas par jour multiple de iphysiq' +! call abort_gcm(modname,abort_message,1) +! ENDIF + + zdtvr = daysec/REAL(day_step) + IF(dtvr.NE.zdtvr) THEN + WRITE(lunout,*) & + 'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr + ENDIF + + ! + ! on remet le calendrier à zero si demande + ! + IF (start_time /= starttime) then + WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' & + ,' fichier restart ne correspond pas à celle lue dans le run.def' + IF (raz_date == 1) then + WRITE(lunout,*)'Je prends l''heure lue dans run.def' + start_time = starttime + ELSE + call abort_gcm("gcm", "'Je m''arrete'", 1) + ENDIF + ENDIF + IF (raz_date == 1) THEN + annee_ref = anneeref + day_ref = dayref + day_ini = dayref + itau_dyn = 0 + itau_phy = 0 + time_0 = 0. + write(lunout,*) & + 'GCM: On reinitialise a la date lue dans gcm.def' + ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN + write(lunout,*) & + 'GCM: Attention les dates initiales lues dans le fichier' + write(lunout,*) & + ' restart ne correspondent pas a celles lues dans ' + write(lunout,*)' gcm.def' + write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref + write(lunout,*)' day_ref=',day_ref," dayref=",dayref + write(lunout,*)' Pas de remise a zero' + ENDIF + +#ifdef CPP_IOIPSL + mois = 1 + heure = 0. + ! Ehouarn: we still need to define JD_ref and JH_ref + ! and since we don't know how many days there are in a year + ! we set JD_ref to 0 (this should be improved ...) + jD_ref=0 + jH_ref=0 +#endif + + if (iflag_phys.eq.1) then + ! these initialisations have already been done (via iniacademic) + ! if running in SW or Newtonian mode + !----------------------------------------------------------------------- + ! Initialisation des constantes dynamiques : + ! ------------------------------------------ + dtvr = zdtvr + CALL iniconst + + !----------------------------------------------------------------------- + ! Initialisation de la geometrie : + ! -------------------------------- + CALL inigeom + + !----------------------------------------------------------------------- + ! Initialisation du filtre : + ! -------------------------- + CALL inifilr + endif ! of if (iflag_phys.eq.1) + ! + !----------------------------------------------------------------------- + ! Initialisation des I/O : + ! ------------------------ + + if (nday>=0) then ! standard case + day_end=day_ini+nday + else ! special case when nday <0, run -nday dynamical steps + day_end=day_ini-nday/day_step + endif + if (less1day) then + day_end=day_ini+floor(time_0+fractday) + endif + if (ndynstep.gt.0) then + day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step)) + endif + + WRITE(lunout,'(a,i7,a,i7)') & + "run from day ",day_ini," to day",day_end + + !----------------------------------------------------------------------- + ! Initialisation de la physique : + ! ------------------------------- + + IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN + ! Physics: +#ifdef CPP_PHYS + CALL iniphysiq(iim,jjm,llm, & + (jjm-1)*iim+2,comm_lmdz, & + daysec,day_ini,dtphys/nsplit_phys, & + rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, & + iflag_phys) +#endif + ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100)) + + + CALL dynredem0("restart.nc", day_end, phis) + ecripar = .TRUE. + +#ifdef CPP_IOIPSL + time_step = zdtvr + if (ok_dyn_ins) then + ! initialize output file for instantaneous outputs + ! t_ops = iecri * daysec ! do operations every t_ops + t_ops =((1.0*iecri)/day_step) * daysec + t_wrt = daysec ! iecri * daysec ! write output every t_wrt + CALL inithist(day_ref,annee_ref,time_step, & + t_ops,t_wrt) + endif + + IF (ok_dyn_ave) THEN + ! initialize output file for averaged outputs + t_ops = iperiod * time_step ! do operations every t_ops + t_wrt = periodav * daysec ! write output every t_wrt + CALL initdynav(day_ref,annee_ref,time_step, & + t_ops,t_wrt) + END IF + dtav = iperiod*dtvr/daysec +#endif +! #endif of #ifdef CPP_IOIPSL + + ! Choix des frequences de stokage pour le offline + ! istdyn=day_step/4 ! stockage toutes les 6h=1jour/4 + ! istdyn=day_step/12 ! stockage toutes les 2h=1jour/12 + istdyn=day_step/4 ! stockage toutes les 6h=1jour/12 + istphy=istdyn/iphysiq + + + ! + !----------------------------------------------------------------------- + ! Integration temporelle du modele : + ! ---------------------------------- + + ! write(78,*) 'ucov',ucov + ! write(78,*) 'vcov',vcov + ! write(78,*) 'teta',teta + ! write(78,*) 'ps',ps + ! write(78,*) 'q',q + + + CALL leapfrog_nogcm_pluto(ucov,vcov,teta,ps,masse,phis,q,time_0) + +END PROGRAM nogcm_pluto + +