! $Id: leapfrog.F90 5105 2024-07-23 17:14:34Z abarral $ ! ! SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0) !IM : pour sortir les param. du modele dans un fis. netcdf 110106 use IOIPSL USE infotrac, ONLY: nqtot, isoCheck USE guide_mod, ONLY: guide_main USE write_field, ONLY: writefield USE control_mod, ONLY: nday, day_step, planet_type, offline, & iconser, iphysiq, iperiod, dissip_period, & iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins, & periodav, ok_dyn_ave, output_grads_dyn use exner_hyb_m, only: exner_hyb use exner_milieu_m, only: exner_milieu USE comvert_mod, ONLY: ap, bp, pressure_exner, presnivs USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf USE logic_mod, ONLY: iflag_phys, ok_guide, forward, leapf, apphys, & statcl, conser, apdiss, purmats, ok_strato USE temps_mod, ONLY: jD_ref, jH_ref, itaufin, day_ini, day_ref, & start_time, dt USE strings_mod, ONLY: msg USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS 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 shema pour l'advection de ! q , en modifiant iadv dans traceur.def (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 "description.h" include "iniprint.h" include "academic.h" REAL, INTENT(IN) :: time_0 ! not used ! 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 :: zqmin, zqmax ! variables dynamiques intermediaire pour le transport REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse ! variables dynamiques au pas -1 REAL :: vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm) REAL :: tetam1(ip1jmp1, llm), psm1(ip1jmp1) REAL :: massem1(ip1jmp1, llm) ! tendances dynamiques REAL :: dv(ip1jm, llm), du(ip1jmp1, llm) REAL :: dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqtot), dp(ip1jmp1) ! tendances de la dissipation REAL :: dvdis(ip1jm, llm), dudis(ip1jmp1, llm) REAL :: dtetadis(ip1jmp1, llm) ! tendances physiques REAL :: dvfi(ip1jm, llm), dufi(ip1jmp1, llm) REAL :: dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqtot), dpfi(ip1jmp1) ! variables pour le fichier histoire REAL :: dtav ! intervalle de temps elementaire REAL :: tppn(iim), tpps(iim), tpn, tps ! INTEGER :: itau, itaufinp1, iav ! INTEGER iday ! jour julien REAL :: time REAL :: SSUM ! REAL finvmaold(ip1jmp1,llm) !ym LOGICAL lafin LOGICAL :: lafin = .FALSE. INTEGER :: ij, iq, l INTEGER :: ik real :: time_step, t_wrt, t_ops ! REAL rdayvrai,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 !IM : pour sortir les param. du modele dans un fis. netcdf 110106 save first data first/.TRUE./ real :: dt_cum character(len = 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 :: physic data callinigrads/.TRUE./ character(len = 10) :: string10 REAL :: flxw(ip1jmp1, llm) ! flux de masse verticale !+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 :: 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(len = 15) :: ztit !IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag. !IM SAVE ip_ebil_dyn !IM DATA ip_ebil_dyn/0/ !-jld character(len = 80) :: dynhist_file, dynhistave_file character(len = *), parameter :: modname = "leapfrog" character(len = 80) :: abort_message logical :: dissip_conservative save dissip_conservative data dissip_conservative/.TRUE./ LOGICAL :: prem save prem DATA prem/.TRUE./ INTEGER :: testita PARAMETER (testita = 9) logical, parameter :: flag_verif = .FALSE. integer :: itau_w ! pas de temps ecriture = itap + itau_phy if (nday>=0) then itaufin = nday * day_step else itaufin = -nday endif itaufinp1 = itaufin + 1 itau = 0 physic = .TRUE. if (iflag_phys==0.or.iflag_phys==2) physic = .FALSE. ! iday = day_ini+itau/day_step ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0 ! IF(time.GT.1.) THEN ! time = time-1. ! iday = iday+1 ! ENDIF !----------------------------------------------------------------------- ! On initialise la pression et la fonction d'Exner : ! -------------------------------------------------- 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 !----------------------------------------------------------------------- ! Debut de l'integration temporelle: ! ---------------------------------- 1 CONTINUE ! Matsuno Forward step begins here ! date: (NB: date remains unchanged for Backward step) ! ----- 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) CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 321') if (ok_guide) then CALL guide_main(itau,ucov,vcov,teta,q,masse,ps) endif ! ! 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 ! ENDIF ! ! 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 ! ... P.Le Van .26/04/94 .... ! Ehouarn: finvmaold is actually not used ! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 ) ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 ) CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 400') 2 CONTINUE ! Matsuno backward or leapfrog step begins here !----------------------------------------------------------------------- ! date: (NB: only leapfrog step requires recomputing date) ! ----- 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 ! gestion des appels de la physique et des dissipations: ! ------------------------------------------------------ ! ! ... P.Le Van ( 6/02/95 ) .... apphys = .FALSE. statcl = .FALSE. conser = .FALSE. apdiss = .FALSE. IF(purmats) THEN ! ! Purely Matsuno time stepping IF(MOD(itau, iconser) ==0.AND. forward) conser = .TRUE. IF(MOD(itau, dissip_period)==0.AND..NOT.forward) & apdiss = .TRUE. IF(MOD(itau, iphysiq)==0.AND..NOT.forward & .and. physic) apphys = .TRUE. ELSE ! ! Leapfrog/Matsuno time stepping IF(MOD(itau, iconser) == 0) conser = .TRUE. IF(MOD(itau + 1, dissip_period)==0 .AND. .NOT. forward) & apdiss = .TRUE. IF(MOD(itau + 1, iphysiq)==0.AND.physic) apphys = .TRUE. END IF ! Ehouarn: for Shallow Water case (ie: 1 vertical layer), ! supress dissipation step if (llm==1) then apdiss = .FALSE. endif CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 589') !----------------------------------------------------------------------- ! calcul des tendances dynamiques: ! -------------------------------- ! ! compute geopotential phi() CALL geopot (ip1jmp1, teta, pk, pks, phis, phi) time = jD_cur + jH_cur CALL caldyn & (itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, & phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time) !----------------------------------------------------------------------- ! calcul des tendances advection des traceurs (dont l'humidite) ! ------------------------------------------------------------- CALL check_isotopes_seq(q, ip1jmp1, & 'leapfrog 686: avant caladvtrac') IF(forward .OR. leapf) THEN ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step CALL caladvtrac(q, pbaru, pbarv, & p, masse, dq, teta, & flxw, pk) ! !write(*,*) 'caladvtrac 346' IF (offline) THEN !maf stokage du flux de masse pour traceurs OFF-LINE CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis, & dtvr, itau) ENDIF ! of IF (offline) ! ENDIF ! of IF( forward .OR. leapf ) !----------------------------------------------------------------------- ! integrations dynamique et traceurs: ! ---------------------------------- CALL msg('720', modname, isoCheck) CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 756') CALL integrd (nqtot, vcovm1, ucovm1, tetam1, psm1, massem1, & dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis) ! $ finvmaold ) CALL msg('724', modname, isoCheck) CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 762') ! .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd) ! !----------------------------------------------------------------------- ! calcul des tendances physiques: ! ------------------------------- ! ######## P.Le Van ( Modif le 6/02/95 ) ########### ! IF(purmats) THEN IF(itau==itaufin.AND..NOT.forward) lafin = .TRUE. ELSE IF(itau + 1 == itaufin) lafin = .TRUE. ENDIF ! ! IF(apphys) THEN ! ! ....... Ajout P.Le Van ( 17/04/96 ) ........... ! 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 ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique ! avec dyn3dmem CALL geopot (ip1jmp1, teta, pk, pks, phis, phi) ! rdaym_ini = itau * dtvr / daysec ! rdayvrai = rdaym_ini + day_ini ! jD_cur = jD_ref + day_ini - day_ref ! $ + int (itau * dtvr / daysec) ! jH_cur = jH_ref + & ! & (itau * dtvr / daysec - int(itau * dtvr / daysec)) jD_cur = jD_ref + day_ini - day_ref + & (itau+1)/day_step IF (planet_type =="generic") THEN ! ! AS: we make jD_cur to be pday jD_cur = int(day_ini + itau / day_step) ENDIF 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) ! 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 ! rajout debug ! lafin = .TRUE. ! Inbterface avec les routines de phylmd (phymars ... ) ! ----------------------------------------------------- !+jld ! Diagnostique de conservation de l'energie : initialisation IF (ip_ebil_dyn>=1) THEN ztit = 'bil dyn' ! Ehouarn: be careful, diagedyn is Earth-specific! IF (planet_type=="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 ) IF (CPPKEY_PHYS) THEN 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) END IF ! ajout des tendances physiques: ! ------------------------------ CALL addfi(dtphys, leapf, forward, & ucov, vcov, teta, q, ps, & dufi, dvfi, dtetafi, dqfi, dpfi) ! ! since addfi updates ps(), also update p(), masse() and pk() 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 IF (ok_strato) THEN CALL top_bound(vcov, ucov, teta, masse, dtphys) ENDIF ! ! Diagnostique de conservation de l'energie : difference IF (ip_ebil_dyn>=1) THEN ztit = 'bil phys' IF (planet_type=="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 ) ENDIF ! of IF( apphys ) IF(iflag_phys==2) THEN ! "Newtonian" case ! Academic case : Simple friction and Newtonan relaxation ! ------------------------------------------------------- DO l = 1, llm DO ij = 1, ip1jmp1 teta(ij, l) = teta(ij, l) - dtvr * & (teta(ij, l) - tetarappel(ij, l)) * (knewt_g + knewt_t(l) * clat4(ij)) ENDDO ENDDO ! of DO l=1,llm if (planet_type=="giant") then ! ! add an intrinsic heat flux at the base of the atmosphere teta(:, 1) = teta(:, 1) + dtvr * aire(:) * ihf / cpp / masse(:, 1) endif CALL friction(ucov, vcov, dtvr) ! ! Sponge layer (if any) IF (ok_strato) THEN ! dufi(:,:)=0. ! dvfi(:,:)=0. ! dtetafi(:,:)=0. ! dqfi(:,:,:)=0. ! dpfi(:)=0. ! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi) CALL top_bound(vcov, ucov, teta, masse, dtvr) ! CALL addfi( dtvr, leapf, forward , ! $ ucov, vcov, teta , q ,ps , ! $ dufi, dvfi, dtetafi , dqfi ,dpfi ) ENDIF ! of IF (ok_strato) ENDIF ! of IF (iflag_phys.EQ.2) !-jld 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 CALL massdair(p, masse) CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1196') !----------------------------------------------------------------------- ! dissipation horizontale et verticale des petites echelles: ! ---------------------------------------------------------- IF(apdiss) THEN ! calcul de l'energie cinetique avant dissipation CALL covcont(llm, ucov, vcov, ucont, vcont) CALL enercin(vcov, ucov, vcont, ucont, ecin0) ! dissipation CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis) ucov = ucov + dudis vcov = vcov + dvdis ! teta=teta+dtetadis !------------------------------------------------------------------------ if (dissip_conservative) then ! On rajoute la tendance due a la transform. Ec -> E therm. cree ! lors de la dissipation CALL covcont(llm, ucov, vcov, ucont, vcont) CALL enercin(vcov, ucov, vcont, ucont, ecin) dtetaecdt = (ecin0 - ecin) / pk ! teta=teta+dtetaecdt dtetadis = dtetadis + dtetaecdt endif teta = teta + dtetadis !------------------------------------------------------------------------ ! ....... P. Le Van ( ajout le 17/04/96 ) ........... ! ... Calcul de la valeur moyenne, unique de h aux poles ..... ! 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 if (1 == 0) then !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics !!! 2) should probably not be here anyway !!! but are kept for those who would want to revert to previous behaviour 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 endif ! of if (1 == 0) END IF ! of IF(apdiss) ! ajout debug ! IF( lafin ) then ! abort_message = 'Simulation finished' ! CALL abort_gcm(modname,abort_message,0) ! ENDIF ! ******************************************************************** ! ******************************************************************** ! .... fin de l'integration dynamique et physique pour le pas itau .. ! ******************************************************************** ! ******************************************************************** ! preparation du pas d'integration suivant ...... CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1509') IF (.NOT.purmats) THEN ! ........................................................ ! .............. schema matsuno + leapfrog .............. ! ........................................................ IF(forward .OR. leapf) THEN itau = itau + 1 ! iday= day_ini+itau/day_step ! time= REAL(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 == 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 !----------------------------------------------------------------------- ! ecriture du fichier histoire moyenne: ! ------------------------------------- IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN IF(itau==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 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q) END IF IF (ok_dyn_ave) THEN CALL writedynav(itau,vcov, & ucov,teta,pk,phi,q,masse,ps,phis) ENDIF ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin)) CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1584') !----------------------------------------------------------------------- ! ecriture de la bande histoire: ! ------------------------------ IF(MOD(itau, iecri)==0) THEN ! ! Ehouarn: output only during LF or Backward Matsuno if (leapf.or.(.not.leapf.and.(.not.forward))) then CALL geopot(ip1jmp1, teta, 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 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) ! 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==itaufin) THEN ! if (planet_type.eq."earth") then ! Write an Earth-format restart file CALL dynredem1("restart.nc", start_time, & vcov, ucov, teta, q, masse, ps) ! endif ! of if (planet_type.eq."earth") CLOSE(99) if (ok_guide) then ! ! set ok_guide to false to avoid extra output ! ! in following forward step ok_guide = .FALSE. endif ! !!! Ehouarn: Why not stop here and now? ENDIF ! of IF (itau.EQ.itaufin) !----------------------------------------------------------------------- ! gestion de l'integration temporelle: ! ------------------------------------ IF(MOD(itau, iperiod)==0) THEN GO TO 1 ELSE IF (MOD(itau - 1, iperiod) == 0) THEN IF(forward) THEN ! fin du pas forward et debut du pas backward forward = .FALSE. leapf = .FALSE. GO TO 2 ELSE ! fin du pas backward et debut du premier pas leapfrog leapf = .TRUE. dt = 2. * dtvr GO TO 2 END IF ! of IF (forward) ELSE ! ...... 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) CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1664') ! ........................................................ ! .............. schema matsuno ............... ! ........................................................ IF(forward) THEN itau = itau + 1 ! iday = day_ini+itau/day_step ! time = REAL(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 == 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 CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1698') IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN IF(itau==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 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q) ENDIF IF (ok_dyn_ave) THEN CALL writedynav(itau,vcov, & ucov,teta,pk,phi,q,masse,ps,phis) ENDIF ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) IF(MOD(itau, iecri)==0) THEN ! IF(MOD(itau,iecri*day_step).EQ.0) THEN CALL geopot(ip1jmp1, teta, 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 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) ! 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==itaufin) THEN ! if (planet_type.eq."earth") then CALL dynredem1("restart.nc", start_time, & vcov, ucov, teta, q, masse, ps) ! endif ! of if (planet_type.eq."earth") if (ok_guide) then ! ! set ok_guide to false to avoid extra output ! ! in following forward step ok_guide = .FALSE. endif ENDIF ! of IF(itau.EQ.itaufin) forward = .TRUE. GO TO 1 ENDIF ! of IF (forward) END IF ! of IF(.not.purmats) END SUBROUTINE leapfrog