--- trunk/libf/dyn3d/leapfrog.f90 2008/03/03 16:32:04 5 +++ trunk/libf/dyn3d/leapfrog.f90 2010/03/03 13:23:49 24 @@ -1,71 +1,65 @@ module leapfrog_m - ! This module is clean: no C preprocessor directive, no include line. - IMPLICIT NONE contains - SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, nq, q, clesphy0, & - time_0) + SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0) ! From dyn3d/leapfrog.F, version 1.6 2005/04/13 08:58:34 - ! Version du 10/01/98, avec coordonnees verticales hybrides, avec - ! nouveaux operat. dissipation * (gradiv2, divgrad2, nxgraro2) + ! Version du 10/01/98, avec coordonnées verticales hybrides, avec + ! nouveaux opérateurs dissipation "*" (gradiv2, divgrad2, nxgraro2) - ! Auteur: P. Le Van /L. Fairhead/F.Hourdin - ! Objet: - ! GCM LMD nouvelle grille + ! Auteurs : P. Le Van, L. Fairhead, F. Hourdin + ! Objet: 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. + ! Dans "inigeom", nouveaux calculs pour les élongations cu, cv + ! et possibilité d'appeler une fonction f(y) à dérivée tangente + ! hyperbolique à la place de la fonction à dérivée sinusoïdale. - ! ... Possibilite de choisir le shema pour l'advection de - ! q, en modifiant iadv dans "traceur.def" (10/02) . + ! Possibilité de choisir le schéma pour l'advection de + ! q, en modifiant iadv dans "traceur.def". - ! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron, 10/99) + ! Pour Van-Leer + vapeur d'eau saturée, iadv(1)=4. ! Pour Van-Leer iadv=10 use dimens_m, only: iim, llm, nqmx - use paramet_m, only: ip1jmp1, ip1jm, llmp1, ijmllm, ijp1llm, jjp1, iip1, & - iip2 + use paramet_m, only: ip1jmp1, ip1jm, ijp1llm, jjp1, iip1 use comconst, only: dtvr, daysec, dtphys use comvert, only: ap, bp use conf_gcm_m, only: day_step, iconser, idissip, iphysiq, iperiod, nday, & offline, periodav - use logic, only: ok_guide, apdiss, apphys, conser, forward, iflag_phys, & - leapf, statcl + use logic, only: ok_guide, iflag_phys use comgeom use serre use temps, only: itaufin, day_ini, dt use iniprint, only: prt_level use com_io_dyn - use abort_gcm_m, only: abort_gcm use ener use calfis_m, only: calfis use exner_hyb_m, only: exner_hyb use guide_m, only: guide use pression_m, only: pression + use pressure_var, only: p3d - integer nq - - INTEGER longcles - PARAMETER (longcles = 20) - REAL clesphy0(longcles) - - ! variables dynamiques + ! Variables dynamiques: REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants REAL teta(ip1jmp1, llm) ! temperature potentielle REAL q(ip1jmp1, llm, nqmx) ! mass fractions of advected fields - REAL ps(ip1jmp1) ! pression au sol - REAL p(ip1jmp1, llmp1) ! pression aux interfac.des couches + REAL ps(ip1jmp1) ! pression au sol, en Pa + REAL masse(ip1jmp1, llm) ! masse d'air + REAL phis(ip1jmp1) ! geopotentiel au sol + + REAL, intent(in):: time_0 + + ! Variables local to the procedure: + + ! Variables dynamiques: + 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) ! geopotential REAL w(ip1jmp1, llm) ! vitesse verticale @@ -93,22 +87,19 @@ REAL tppn(iim), tpps(iim), tpn, tps - INTEGER itau, itaufinp1 + INTEGER itau ! index of the time step of the dynamics, starts at 0 INTEGER iday ! jour julien - REAL time ! Heure de la journee en fraction d'1 jour + REAL time ! time of day, as a fraction of day length REAL SSUM - REAL time_0, finvmaold(ip1jmp1, llm) + real finvmaold(ip1jmp1, llm) LOGICAL :: lafin=.false. INTEGER ij, l REAL rdayvrai, rdaym_ini - LOGICAL callinigrads - data callinigrads/.true./ - - !+jld variables test conservation energie + ! 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 @@ -116,28 +107,16 @@ REAL dtetaecdt(ip1jmp1, llm) REAL vcont(ip1jm, llm), ucont(ip1jmp1, llm) CHARACTER*15 ztit - INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag. - SAVE ip_ebil_dyn - DATA ip_ebil_dyn /0/ - - character(len=*), parameter:: modname = "leapfrog" - character*80 abort_message - - logical dissip_conservative - save dissip_conservative - data dissip_conservative /.true./ - - LOGICAL prem - save prem - DATA prem /.true./ + INTEGER:: ip_ebil_dyn = 0 ! PRINT level for energy conserv. diag. + + logical:: dissip_conservative = .true. + logical forward, leapf, apphys, conser, apdiss !--------------------------------------------------- print *, "Call sequence information: leapfrog" itaufin = nday * day_step - itaufinp1 = itaufin + 1 - itau = 0 iday = day_ini time = time_0 @@ -146,21 +125,21 @@ iday = iday + 1 ENDIF + dq = 0. ! On initialise la pression et la fonction d'Exner : - dq=0. - CALL pression(ip1jmp1, ap, bp, ps, p) - CALL exner_hyb(ps, p, pks, pk, pkf) + CALL pression(ip1jmp1, ap, bp, ps, p3d) + CALL exner_hyb(ps, p3d, pks, pk, pkf) ! Debut de l'integration temporelle: - do - if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600) then + outer_loop:do + if (ok_guide.and.(itaufin - itau - 1) * dtvr > 21600.) then call guide(itau, ucov, vcov, teta, q, masse, ps) else IF (prt_level > 9) print *, & - 'Attention : on ne guide pas les 6 dernieres heures.' + 'Attention : on ne guide pas les 6 dernières heures.' endif - CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1) + CALL SCOPY(ip1jm * llm, 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) @@ -177,7 +156,6 @@ ! gestion des appels de la physique et des dissipations: apphys = .FALSE. - statcl = .FALSE. conser = .FALSE. apdiss = .FALSE. @@ -196,7 +174,7 @@ ! calcul des tendances advection des traceurs (dont l'humidite) IF (forward .OR. leapf) THEN - CALL caladvtrac(q, pbaru, pbarv, p, masse, dq, teta, pk) + CALL caladvtrac(q, pbaru, pbarv, p3d, masse, dq, teta, pk) IF (offline) THEN !maf stokage du flux de masse pour traceurs OFF-LINE CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, dtvr, & @@ -206,15 +184,16 @@ ! integrations dynamique et traceurs: CALL integrd(2, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, & - dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, finvmaold) + dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis, & + finvmaold, leapf) ! calcul des tendances physiques: IF (apphys) THEN IF (itau + 1 == itaufin) lafin = .TRUE. - CALL pression(ip1jmp1, ap, bp, ps, p) - CALL exner_hyb(ps, p, pks, pk, pkf) + CALL pression(ip1jmp1, ap, bp, ps, p3d) + CALL exner_hyb(ps, p3d, pks, pk, pkf) rdaym_ini = itau * dtvr / daysec rdayvrai = rdaym_ini + day_ini @@ -224,13 +203,13 @@ ! Diagnostique de conservation de l'énergie : initialisation IF (ip_ebil_dyn >= 1) THEN ztit='bil dyn' - CALL diagedyn(ztit, 2, 1, 1, dtphys & - , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2)) + CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, & + teta, q(:, :, 1), q(:, :, 2)) ENDIF - CALL calfis(nq, lafin, rdayvrai, time, ucov, vcov, teta, q, & - masse, ps, p, pk, phis, phi, du, dv, dteta, dq, w, & - clesphy0, dufi, dvfi, dtetafi, dqfi, dpfi) + CALL calfis(nqmx, lafin, rdayvrai, time, ucov, vcov, teta, q, & + masse, ps, pk, phis, phi, du, dv, dteta, dq, w, & + dufi, dvfi, dtetafi, dqfi, dpfi) ! ajout des tendances physiques: CALL addfi(nqmx, dtphys, & @@ -240,13 +219,13 @@ ! Diagnostique de conservation de l'énergie : difference IF (ip_ebil_dyn >= 1) THEN ztit = 'bil phys' - CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p, pk, & + CALL diagedyn(ztit, 2, 1, 1, dtphys, ucov, vcov, ps, p3d, pk, & teta, q(:, :, 1), q(:, :, 2)) ENDIF ENDIF - CALL pression(ip1jmp1, ap, bp, ps, p) - CALL exner_hyb(ps, p, pks, pk, pkf) + CALL pression(ip1jmp1, ap, bp, ps, p3d) + CALL exner_hyb(ps, p3d, pks, pk, pkf) ! dissipation horizontale et verticale des petites echelles: @@ -256,7 +235,7 @@ call enercin(vcov, ucov, vcont, ucont, ecin0) ! dissipation - CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis) + CALL dissip(vcov, ucov, teta, p3d, dvdis, dudis, dtetadis) ucov=ucov + dudis vcov=vcov + dvdis @@ -315,10 +294,7 @@ ENDIF ENDIF - IF (itau == itaufinp1) then - abort_message = 'Simulation finished' - call abort_gcm(modname, abort_message, 0) - ENDIF + IF (itau == itaufin + 1) exit outer_loop ! ecriture du fichier histoire moyenne: @@ -332,7 +308,7 @@ ENDIF IF (itau == itaufin) THEN - CALL dynredem1("restart.nc", 0., vcov, ucov, teta, q, masse, ps) + CALL dynredem1("restart.nc", vcov, ucov, teta, q, masse, ps) CLOSE(99) ENDIF @@ -350,12 +326,12 @@ dt = 2. * dtvr END IF ELSE - ! ...... pas leapfrog ..... + ! pas leapfrog leapf = .TRUE. dt = 2. * dtvr END IF end do - end do + end do outer_loop END SUBROUTINE leapfrog