! ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/physiq.F,v 1.8 2005/02/24 09:58:18 fairhead Exp $ ! c SUBROUTINE physiq (nlon,nlev,nqmax, . debut,lafin,rjourvrai,gmtime,pdtphys, . paprs,pplay,ppk,pphi,pphis,presnivs, . u,v,t,qx, . omega, . d_u, d_v, d_t, d_qx, d_ps) c====================================================================== c c Modifications pour la physique de Titan c adaptation a partir de celle de Venus c S. Lebonnois (LMD/CNRS) Mai 2008 c c --------------------------------------------------------------------- c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818 c c Objet: Moniteur general de la physique du modele cAA Modifications quant aux traceurs : cAA - uniformisation des parametrisations ds phytrac cAA - stockage des moyennes des champs necessaires cAA en mode traceur off-line c modif ( P. Le Van , 12/10/98 ) c c Arguments: c c nlon------input-I-nombre de points horizontaux c nlev------input-I-nombre de couches verticales c nqmax-----input-I-nombre de traceurs c debut-----input-L-variable logique indiquant le premier passage c lafin-----input-L-variable logique indiquant le dernier passage c rjourvrai-input-R-NBjours c gmtime----input-R-temps universel dans la journee (fraction de jour) c pdtphys---input-R-pas d'integration pour la physique (seconde) c paprs-----input-R-pression pour chaque inter-couche (en Pa) c pplay-----input-R-pression pour le mileu de chaque couche (en Pa) c ppk -----input-R-fonction d'Exner au milieu de couche c pphi------input-R-geopotentiel de chaque couche (g z) (reference sol) c pphis-----input-R-geopotentiel du sol c presnivs--input_R_pressions approximat. des milieux couches ( en PA) c u---------input-R-vitesse dans la direction X (de O a E) en m/s c v---------input-R-vitesse Y (de S a N) en m/s c t---------input-R-temperature (K) c qx--------input-R-mass mixing ratio traceurs (kg/kg) c d_t_dyn---input-R-tendance dynamique pour "t" (K/s) c omega-----input-R-vitesse verticale en Pa/s c c d_u-----output-R-tendance physique de "u" (m/s/s) c d_v-----output-R-tendance physique de "v" (m/s/s) c d_t-----output-R-tendance physique de "t" (K/s) c d_qx----output-R-tendance physique de "qx" (kg/kg/s) c d_ps----output-R-tendance physique de la pression au sol c====================================================================== USE ioipsl USE histcom USE infotrac USE control_mod use dimphy USE comgeomphy IMPLICIT none c====================================================================== c CLEFS CPP POUR LES IO c ===================== #define histmth #define histday #define histins c====================================================================== #include "dimensions.h" integer jjmp1 parameter (jjmp1=jjm+1-1/jjm) #include "dimsoil.h" #include "clesphys.h" #include "temps.h" #include "iniprint.h" #include "timerad.h" #include "logic.h" #include "comorbit.h" #include "microtab.h" #include "diagmuphy.h" #include "itemps.h" c====================================================================== LOGICAL ok_mensuel ! sortir le fichier mensuel save ok_mensuel c PARAMETER (ok_mensuel=.true.) c LOGICAL ok_journe ! sortir le fichier journalier save ok_journe c PARAMETER (ok_journe=.true.) c LOGICAL ok_instan ! sortir le fichier instantane save ok_instan c PARAMETER (ok_instan=.true.) c c====================================================================== c c Variables argument: c INTEGER nlon INTEGER nlev INTEGER nqmax REAL rjourvrai REAL gmtime REAL pdtphys LOGICAL debut, lafin REAL paprs(klon,klev+1) REAL pplay(klon,klev) REAL pphi(klon,klev) REAL pphis(klon) REAL presnivs(klev) ! ADAPTATION GCM POUR CP(T) REAL ppk(klon,klev) REAL u(klon,klev) REAL v(klon,klev) REAL t(klon,klev) REAL qx(klon,klev,nqmax) REAL,save,allocatable :: t_ancien(:,:) REAL,save,allocatable :: u_ancien(:,:) LOGICAL ancien_ok SAVE ancien_ok REAL d_u_dyn(klon,klev) REAL d_t_dyn(klon,klev) REAL omega(klon,klev) REAL d_u(klon,klev) REAL d_v(klon,klev) REAL d_t(klon,klev) REAL d_qx(klon,klev,nqmax) REAL d_ps(klon) REAL,save,allocatable :: swnet(:,:) REAL,save,allocatable :: lwnet(:,:) c c Variables propres a la physique c REAL dtime SAVE dtime ! pas temporel de la physique c INTEGER radpas SAVE radpas ! frequence d'appel rayonnement c INTEGER chimpas SAVE chimpas ! frequence d'appel chimie c REAL,save,allocatable :: radsol(:) ! bilan radiatif au sol calcule par code radiatif REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche) INTEGER,save :: itap ! compteur pour la physique REAL,save,allocatable :: ftsol(:) ! temperature du sol REAL,save,allocatable :: ftsoil(:,:) ! temperature dans le sol REAL,save,allocatable :: falbe(:) ! albedo REAL delp(klon,klev) ! epaisseur d'une couche CMODDEB FLOTT c c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM): c REAL,save,allocatable :: zmea(:) ! orographie moyenne REAL,save,allocatable :: zstd(:) ! deviation standard de l'OESM REAL,save,allocatable :: zsig(:) ! pente de l'OESM REAL,save,allocatable :: zgam(:) ! anisotropie de l'OESM REAL,save,allocatable :: zthe(:) ! orientation de l'OESM REAL,save,allocatable :: zpic(:) ! Maximum de l'OESM REAL,save,allocatable :: zval(:) ! Minimum de l'OESM REAL,save,allocatable :: rugoro(:) ! longueur de rugosite de l'OESM INTEGER igwd,idx(klon),itest(klon) c c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro REAL zulow(klon),zvlow(klon) REAL zustrdr(klon), zvstrdr(klon) REAL zustrli(klon), zvstrli(klon) REAL zustrhi(klon), zvstrhi(klon) c Pour calcul GW drag oro et nonoro: CALCUL de N2: real zdzlev(klon,klev) real ztlev(klon,klev),zpklev(klon,klev) real ztetalay(klon,klev),ztetalev(klon,klev) real zdtetalev(klon,klev) real zn2(klon,klev) ! BV^2 at plev c Pour les bilans de moment angulaire, integer bilansmc c Pour le transport de ballons integer ballons c j'ai aussi besoin c du stress de couche limite a la surface: REAL zustrcl(klon),zvstrcl(klon) c et du stress total c de la physique: REAL zustrph(klon),zvstrph(klon) c REAL,save,allocatable :: zuthe(:),zvthe(:) c Variables locales: c REAL cdragh(klon) ! drag coefficient pour T and Q REAL cdragm(klon) ! drag coefficient pour vent c cAA Pour TRACEURS cAA REAL,save,allocatable :: source(:,:) integer nmicro save nmicro REAL ycoefh(klon,klev) ! coef d'echange pour phytrac REAL yu1(klon) ! vents dans la premiere couche U REAL yv1(klon) ! vents dans la premiere couche V character*8 nom REAL qaer(klon,klev,nqmax) REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee REAL,save,allocatable :: dlw(:) ! derivee infra rouge REAL,save,allocatable :: fder(:) ! Derive de flux (sensible et latente) REAL ve(klon) ! integr. verticale du transport meri. de l'energie REAL vq(klon) ! integr. verticale du transport meri. de l'eau REAL ue(klon) ! integr. verticale du transport zonal de l'energie REAL uq(klon) ! integr. verticale du transport zonal de l'eau c c====================================================================== c c Declaration des procedures appelees c EXTERNAL ajsec ! ajustement sec EXTERNAL clmain ! couche limite EXTERNAL hgardfou ! verifier les temperatures EXTERNAL orbite ! calculer l'orbite EXTERNAL phyetat0 ! lire l'etat initial de la physique EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique EXTERNAL radlwsw ! rayonnements solaire et infrarouge EXTERNAL suphec ! initialiser certaines constantes c EXTERNAL transp ! transport total de l'eau et de l'energie EXTERNAL abort_gcm EXTERNAL printflag EXTERNAL zenang EXTERNAL diagetpq EXTERNAL conf_phys EXTERNAL diagphy EXTERNAL mucorr EXTERNAL phytrac c c Variables locales c CXXX PB REAL fluxt(klon,klev) ! flux turbulent de chaleur REAL fluxu(klon,klev) ! flux turbulent de vitesse u REAL fluxv(klon,klev) ! flux turbulent de vitesse v c REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec REAL flux_ec(klon,klev) ! flux de chaleur Ec c c Le rayonnement n'est pas calcule tous les pas, il faut donc c sauvegarder les sorties du rayonnement REAL,save,allocatable :: heat(:,:) ! chauffage solaire REAL,save,allocatable :: cool(:,:) ! refroidissement infrarouge REAL,save,allocatable :: dtrad(:,:) ! K s-1 REAL,save,allocatable :: topsw(:), toplw(:) REAL,save,allocatable :: solsw(:), sollw(:) REAL,save,allocatable :: sollwdown(:) ! downward LW flux at surface REAL tmpout(klon,klev) ! K s-1 INTEGER itaprad SAVE itaprad REAL zdtime c c CHIMIE REAL dtimechim INTEGER itapchim,appel_chim SAVE itapchim,dtimechim c ORBITE REAL dist, rmu0(klon), fract(klon), pdecli REAL zday REAL zls c INTEGER i, k, iq, ig, j, ll, l c REAL zphi(klon,klev) REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2 c c Variables du changement c c ajs: ajustement sec c vdf: couche limite (Vertical DiFfusion) c mph: microphysique c kim: chimie REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax) REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev) c REAL d_ts(klon) c REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev) REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax) c REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax) CMOD LOTT: Tendances Orography Sous-maille REAL d_u_oro(klon,klev), d_v_oro(klon,klev) REAL d_t_oro(klon,klev) REAL d_u_lif(klon,klev), d_v_lif(klon,klev) REAL d_t_lif(klon,klev) C Tendances Ondes de G non oro (runs strato). REAL d_u_hin(klon,klev), d_v_hin(klon,klev) REAL d_t_hin(klon,klev) c c Variables liees a l'ecriture de la bande histoire physique c INTEGER ecrit_mth SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel) c INTEGER ecrit_day SAVE ecrit_day ! frequence d'ecriture (fichier journalier) c INTEGER ecrit_ins SAVE ecrit_ins ! frequence d'ecriture (fichier instantane) c integer itau_w ! pas de temps ecriture = itap + itau_phy c Variables locales pour effectuer les appels en serie c REAL t_seri(klon,klev) REAL u_seri(klon,klev), v_seri(klon,klev) c REAL tr_seri(klon,klev,nqmax) INTEGER length PARAMETER ( length = 100 ) REAL tabcntr0( length ) c c pour ioipsl INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev) REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D REAL zx_tmp_2d(iim,jjmp1),zx_tmp_3d(iim,jjmp1,klev) REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1) INTEGER nid_day, nid_mth, nid_ins SAVE nid_day, nid_mth, nid_ins c INTEGER nhori, nvert, idayref REAL zsto, zout, zsto1, zsto2, zero parameter (zero=0.0e0) real zjulian save zjulian CHARACTER*2 str2 character*20 modname character*80 abort_message logical ok_sync character*30 nom_fichier character*10 varname character*40 vartitle character*20 varunits C Variables liees au bilan d'energie et d'enthalpi REAL ztsol(klon) REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot $ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot $ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec REAL d_h_vcol_phy REAL fs_bound, fq_bound SAVE d_h_vcol_phy REAL zero_v(klon),zero_v2(klon,klev) CHARACTER*15 ztit INTEGER ip_ebil ! PRINT level for energy conserv. diag. SAVE ip_ebil DATA ip_ebil/2/ INTEGER if_ebil ! level for energy conserv. dignostics SAVE if_ebil c+jld ec_conser REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique c-jld ec_conser c TEST VENUS... REAL mang(klon,klev) ! moment cinetique REAL mangtot ! moment cinetique total CHARACTER*2 str1 c Temporaire avant de trouver mieux : c Recuperation des TAU du TR REAL t_tauhvd(klon,klev),t_khvd(klon,klev) REAL t_tcld(klon,klev),t_kcld(klon,klev) REAL t_kcvd(klon,klev) c ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX INTEGER ngrid PARAMETER (ngrid=(jjm-1)*iim+2) ! = klon INTEGER NSPECV,NSPECI,NLAYER PARAMETER (NSPECV=24,NSPECI=46,NLAYER=llm) REAL TAUHID(ngrid,NLAYER,NSPECI) & ,TAUCID(ngrid,NLAYER,NSPECI) & ,TAUGID(ngrid,NLAYER,NSPECI) & ,TAUHVD(ngrid,NLAYER,NSPECV) & ,TAUCVD(ngrid,NLAYER,NSPECV) & ,TAUGVD(ngrid,NLAYER,NSPECV) COMMON /TAUD/ TAUHID,TAUCID,TAUGID,TAUHVD,TAUCVD,TAUGVD * common relatifs au nuages real rmcbar(ngrid,NLAYER),xfbar(ngrid,NLAYER,4) integer ncount(ngrid,NLAYER) COMMON/rnuabar/ncount,rmcbar,xfbar REAL ch4(klon,jjm+1),dch4(jjm+1) INTEGER ig0 integer ich4 common/ch4ind/ich4 c flux de chaleur latente d'evaporation CH4 REAL fclat(klon) c reservoir de surface REAL,save,allocatable :: reservoir(:) c Declaration des constantes et des fonctions thermodynamiques c #include "YOMCST.h" c====================================================================== c INITIALISATIONS c====================================================================== modname = 'physiq' ok_sync=.TRUE. bilansmc = 0 ballons = 0 IF (if_ebil.ge.1) THEN DO i=1,klon zero_v(i)=0. END DO DO i=1,klon DO j=1,klev zero_v2(i,j)=0. END DO END DO END IF c PREMIER APPEL SEULEMENT c======================== IF (debut) THEN allocate(t_ancien(klon,klev),u_ancien(klon,klev)) allocate(swnet(klon,klevp1),lwnet(klon,klevp1)) allocate(radsol(klon),ftsol(klon),falbe(klon)) allocate(rlev(klon,klevp1),ftsoil(klon,nsoilmx)) allocate(zmea(klon),zstd(klon),zsig(klon),zgam(klon)) allocate(zthe(klon),zpic(klon),zval(klon),rugoro(klon)) allocate(zuthe(klon),zvthe(klon),dlw(klon),fder(klon)) allocate(heat(klon,klev),cool(klon,klev)) allocate(dtrad(klon,klev),topsw(klon),toplw(klon)) allocate(solsw(klon),sollw(klon),sollwdown(klon)) allocate(source(klon,nqmax)) allocate(reservoir(klon)) CALL suphec ! initialiser constantes et parametres phys. IF (if_ebil.ge.1) d_h_vcol_phy=0. c c appel a la lecture du physiq.def c call conf_phys(ok_mensuel,ok_journe,ok_instan,if_ebil) c c Initialiser les compteurs: c itap = 0 itaprad = 0 itapchim = 0 ncount(:,:) = 0 c c Lecture startphy.nc : c c REMETTRE TOUS LES PARAMETRES POUR OROGW... A FAIRE POUR TITAN CALL phyetat0 ("startphy.nc",dtime, . rlatd,rlond,ftsol,ftsoil, . falbe, solsw, sollw, . dlw,radsol,reservoir, c . zmea,zstd,zsig,zgam,zthe,zpic,zval, . tabcntr0, . t_ancien, ancien_ok) c dtime est lu dans startphy c pdtphys est calcule a partir des nouvelles conditions: c Reinitialisation du pas de temps physique quand changement IF (ABS(dtime-pdtphys).GT.0.001) THEN WRITE(lunout,*) 'Pas physique a change',dtime, . pdtphys c abort_message='Pas physique n est pas correct ' c call abort_gcm(modname,abort_message,1) c---------------- c pour initialiser convenablement le time_counter, il faut tenir compte c du changement de dtime en changeant itau_phy (point de depart) itau_phy = NINT(itau_phy*dtime/pdtphys) c---------------- dtime=pdtphys ENDIF radpas = NINT( RDAY/pdtphys/nbapp_rad) chimpas = radpas*nbapp_rad/nbapp_chim CALL printflag( tabcntr0,radpas,chimpas,ok_mensuel, . ok_journe, ok_instan ) c c Initialiser les pas de temps: c dtimerad = dtime*FLOAT(radpas) ! pas de temps du rayonnement (s) c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas dtimechim = dtime*FLOAT(chimpas) ! pas de temps de la chimie (s) c PRINT*,'dtimechim,dtime,chimpas',dtimechim,dtime,chimpas c INITIALISATION ORBITE CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) c--------- c FLOTT c IF (ok_orodr) THEN c DO i=1,klon c rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0) c ENDDO c CALL SUGWD(klon,klev,paprs,pplay) c DO i=1,klon c zuthe(i)=0. c zvthe(i)=0. c if(zstd(i).gt.10.)then c zuthe(i)=(1.-zgam(i))*cos(zthe(i)) c zvthe(i)=(1.-zgam(i))*sin(zthe(i)) c endif c ENDDO c ENDIF if (bilansmc.eq.1) then C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES C DU BILAN DE MOMENT ANGULAIRE. open(27,file='aaam_bud.out',form='formatted') open(28,file='fields_2d.out',form='formatted') write(*,*)'Ouverture de aaam_bud.out (FL Vous parle)' write(*,*)'Ouverture de fields_2d.out (FL Vous parle)' endif !bilansmc c--------------SLEBONNOIS C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES C DES BALLONS if (ballons.eq.1) then open(30,file='ballons-lat.out',form='formatted') open(31,file='ballons-lon.out',form='formatted') open(32,file='ballons-u.out',form='formatted') open(33,file='ballons-v.out',form='formatted') open(34,file='ballons-alt.out',form='formatted') write(*,*)'Ouverture des ballons*.out' endif !ballons c------------- c--------- C TRACEURS C source dans couche limite source = 0.0 ! pas de source, pour l'instant C c Si microphysique offline, pas besoin d'avoir de traceurs microphysiques c car on lit les profils verticaux des qaer dans une look-up table pour c le rayonnement. c calcul de nmicro c !!!! Les traceurs microphysiques doivent etre toujours en premiers!! nmicro = 0 do iq=1,nqmax nom = tname(iq) c print*,iq,"nom=",nom,"tname=",tname(iq) print*,iq,"nom=",nom if (nom(1:1).eq."q") then nmicro = nmicro+1 endif enddo print*,"nmicro=",nmicro c c Verifications: c IF ((nmicro.eq.0).and.(microfi.eq.1)) THEN abort_message="MICROPHYSIQUE ONLINE, MAIS NMICRO=0..." call abort_gcm(modname,abort_message,1) ENDIF IF (microfi.lt.1.and.clouds.eq.1) THEN write(lunout,*)"microfi.lt.1.and.clouds.eq.1" abort_message = & "Impossible de faire des nuages sans microphysique..." call abort_gcm(modname,abort_message,1) ENDIF IF (nlon .NE. klon) THEN WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, . klon abort_message='nlon et klon ne sont pas coherents' call abort_gcm(modname,abort_message,1) ENDIF IF (nlev .NE. klev) THEN WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, . klev abort_message='nlev et klev ne sont pas coherents' call abort_gcm(modname,abort_message,1) ENDIF c IF (dtime*FLOAT(radpas).GT.(RDAY*0.25).AND.cycle_diurne) $ THEN WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant' WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne" abort_message='Nbre d appels au rayonnement insuffisant' call abort_gcm(modname,abort_message,1) ENDIF c WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=", . iflag_ajs c ecrit_mth = NINT(RDAY/dtime) *nday ! tous les nday jours IF (ok_mensuel) THEN WRITE(lunout,*)'La frequence de sortie mensuelle est de ', . ecrit_mth ENDIF ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours IF (ok_journe) THEN WRITE(lunout,*)'La frequence de sortie journaliere est de ', . ecrit_day ENDIF ecrit_ins = NINT(RDAY/dtime*ecritphy) ! Fraction de jour reglable IF (ok_instan) THEN WRITE(lunout,*)'La frequence de sortie instant. est de ', . ecrit_ins ENDIF c Initialisation des sorties c=========================== #ifdef CPP_IOIPSL #ifdef histmth #include "ini_histmth.h" #endif #ifdef histday #include "ini_histday.h" #endif #ifdef histins #include "ini_histins.h" #endif #endif c c Initialiser les valeurs de u pour calculs tendances c (pour T, c'est fait dans phyetat0) c DO k = 1, klev DO i = 1, klon u_ancien(i,k) = u(i,k) ENDDO ENDDO rmcbar = 0. xfbar = 0. ENDIF ! debut c==================================================================== c====================================================================== c Creer un reservoir de surface infini c reservoir(:) = 2. c Mettre a zero des variables de sortie (pour securite) c DO i = 1, klon d_ps(i) = 0.0 ENDDO DO k = 1, klev DO i = 1, klon d_t(i,k) = 0.0 d_u(i,k) = 0.0 d_v(i,k) = 0.0 ENDDO ENDDO DO iq = 1, nqmax DO k = 1, klev DO i = 1, klon d_qx(i,k,iq) = 0.0 ENDDO ENDDO ENDDO c c Ne pas affecter les valeurs entrees de u, v, h, et q c DO k = 1, klev DO i = 1, klon t_seri(i,k) = t(i,k) u_seri(i,k) = u(i,k) v_seri(i,k) = v(i,k) ENDDO ENDDO DO iq = 1, nqmax DO k = 1, klev DO i = 1, klon tr_seri(i,k,iq) = qx(i,k,iq) ENDDO ENDDO ENDDO C DO i = 1, klon ztsol(i) = ftsol(i) ENDDO C IF (if_ebil.ge.1) THEN ztit='after dynamic' CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) C Comme les tendances de la physique sont ajoute dans la dynamique, C on devrait avoir que la variation d'entalpie par la dynamique C est egale a la variation de la physique au pas de temps precedent. C Donc la somme de ces 2 variations devrait etre nulle. call diagphy(airephy,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, zero_v e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol+d_h_vcol_phy, d_qt, 0. s , fs_bound, fq_bound ) END IF c==================================================================== c Diagnostiquer la tendance dynamique c IF (ancien_ok) THEN DO k = 1, klev DO i = 1, klon d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime ENDDO ENDDO ! ADAPTATION GCM POUR CP(T) do i=1,klon flux_dyn(i,1) = 0.0 do j=2,klev flux_dyn(i,j) = flux_dyn(i,j-1) . +cpdet(t_seri(i,j-1))/RG*d_t_dyn(i,j-1)*(paprs(i,j-1)-paprs(i,j)) enddo enddo ELSE DO k = 1, klev DO i = 1, klon d_u_dyn(i,k) = 0.0 d_t_dyn(i,k) = 0.0 ENDDO ENDDO ancien_ok = .TRUE. ENDIF c==================================================================== c c Ajouter le geopotentiel du sol: c DO k = 1, klev DO i = 1, klon zphi(i,k) = pphi(i,k) + pphis(i) ENDDO ENDDO c calcul du geopotentiel aux niveaux intercouches c ponderation des altitudes au niveau des couches en dp/p DO l=1,klev DO i=1,klon c zzlay(i,l)=zphi(i,l)/RG c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE: zzlay(i,l)=RG*RA*RA/(RG*RA-zphi(i,l))-RA ENDDO ENDDO DO i=1,klon c zzlev(i,1)=0. c CORRECTION 13/01/2011 c (correspond a la position de la surface en ce point vs RA) zzlev(i,1)=pphis(i)/RG ENDDO DO l=2,klev DO i=1,klon z1=(pplay(i,l-1)+paprs(i,l))/(pplay(i,l-1)-paprs(i,l)) z2=(paprs(i,l)+pplay(i,l))/(paprs(i,l)-pplay(i,l)) zzlev(i,l)=(z1*zzlay(i,l-1)+z2*zzlay(i,l))/(z1+z2) ENDDO ENDDO DO i=1,klon zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev)) ENDDO c- - - - - - - - - - - - - - - - c DIAGNOSTIQUE GRILLE VERTICALE c- - - - - - - - - - - - - - - - c print*,"DIAGNOSTIQUE GRILLE VERTICALE" c i=klon/2 c print*,"Niveau Pression Altitude (lev puis lay)" c do l=1,klev c print*,l,paprs(i,l),zzlev(i,l) c print*,l,pplay(i,l),zzlay(i,l) c enddo c print*,klev+1,paprs(i,klev+1),zzlev(i,klev+1) c stop c==================================================================== c c Verifier les temperatures c CALL hgardfou(t_seri,ftsol,'debutphy') c==================================================================== c c Incrementer le compteur de la physique c itap = itap + 1 c==================================================================== c c Epaisseurs couches DO k = 1, klev DO i = 1, klon delp(i,k) = paprs(i,k)-paprs(i,k+1) ENDDO ENDDO c==================================================================== c ORBITE ET ECLAIREMENT c==================================================================== c Pour TITAN: c calcul de la longitude solaire CALL solarlong(rjourvrai+gmtime,zls) print*,'Ls',zls*180./RPI ! zls est en radians !! CALL orbite(zls,dist,pdecli) c dans zenang, Ls en degres ; dans mucorr, Ls en radians IF (cycle_diurne) THEN zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s) CALL zenang(zls*180./RPI,gmtime,zdtime,rlatd,rlond,rmu0,fract) ELSE call mucorr(klon,zls,rlatd,rmu0,fract) ENDIF c==================================================================== c COUCHE LIMITE c==================================================================== c------------------------------- c TEST: on ne tient pas compte des calculs de clmain mais on force c l'equilibre radiatif du sol if (1.eq.0) then if (debut) then print*,"ATTENTION, CLMAIN SHUNTEE..." endif DO i = 1, klon sens(i) = 0.0e0 ! flux de chaleur sensible au sol fder(i) = 0.0e0 dlw(i) = 0.0e0 ENDDO c Incrementer la temperature du sol c DO i = 1, klon d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200 ftsol(i) = ftsol(i) + d_ts(i) do j=1,nsoilmx ftsoil(i,j)=ftsol(i) enddo ENDDO c------------------------------- else c------------------------------- fder = dlw c print*,"radsol avant clmain=",radsol(klon/2) c print*,"solsw avant clmain=",solsw(klon/2) c print*,"sollw avant clmain=",sollw(klon/2) c CLMAIN ! ADAPTATION GCM POUR CP(T) CALL clmain(dtime,itap, e t_seri,u_seri,v_seri, e rmu0, e ftsol, $ ftsoil, $ paprs,pplay,ppk,radsol,falbe, e solsw, sollw, sollwdown, fder, e rlond, rlatd, cuphy, cvphy, e debut, lafin, s d_t_vdf,d_u_vdf,d_v_vdf,d_ts, s fluxt,fluxu,fluxv,cdragh,cdragm, s dsens, s ycoefh,yu1,yv1) c print*,"radsol apres clmain=",radsol(klon/2) c print*,"solsw apres clmain=",solsw(klon/2) c print*,"sollw apres clmain=",sollw(klon/2) CXXX Incrementation des flux DO i = 1, klon sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol fder(i) = dlw(i) + dsens(i) ENDDO CXXX DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k) d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k) d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k) d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s ENDDO ENDDO c print*,"d_t_vdf1=",d_t_vdf(1,:)*dtime c print*,"d_t_vdf2=",d_t_vdf(klon/2,:)*dtime c print*,"d_t_vdf3=",d_t_vdf(klon,:)*dtime c print*,"d_u_vdf=",d_u_vdf(klon/2,:)*dtime c print*,"d_v_vdf=",d_v_vdf(klon/2,:)*dtime C TRACEURS d_tr_vdf = 0. if (iflag_trac.eq.1) then DO iq=1, nqmax CALL cltrac(dtime,ycoefh,t_seri, s tr_seri(1,1,iq), source, e paprs, pplay, delp, s d_tr_vdf(1,1,iq)) tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq) d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s ENDDO endif IF (if_ebil.ge.2) THEN ztit='after clmain' CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(airephy,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, sens e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) END IF C c c Incrementer la temperature du sol c c print*,'Tsol avant clmain:',ftsol(klon/2) DO i = 1, klon ftsol(i) = ftsol(i) + d_ts(i) ENDDO c print*,'DTsol apres clmain:',d_ts(klon/2) c print*,'Tsol apres clmain:',ftsol(klon/2) c Calculer la derive du flux infrarouge c DO i = 1, klon dlw(i) = - 4.0*RSIGMA*ftsol(i)**3 ENDDO c------------------------------- endif ! fin du TEST c c Appeler l'ajustement sec c c=================================================================== c Convection seche c=================================================================== c d_t_ajs(:,:)=0. d_u_ajs(:,:)=0. d_v_ajs(:,:)=0. d_tr_ajs(:,:,:)=0. c IF(prt_level>9)WRITE(lunout,*) . 'AVANT LA CONVECTION SECHE , iflag_ajs=' s ,iflag_ajs if(iflag_ajs.eq.0) then c Rien c ==== IF(prt_level>9)WRITE(lunout,*)'pas de convection' else if(iflag_ajs.eq.1) then c Ajustement sec c ============== IF(prt_level>9)WRITE(lunout,*)'ajsec' ! ADAPTATION GCM POUR CP(T) CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax, . tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs) ! ADAPTATION GCM POUR CP(T) do i=1,klon flux_ajs(i,1) = 0.0 do j=2,klev flux_ajs(i,j) = flux_ajs(i,j-1) . + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime . *delp(i,j-1) enddo enddo t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:) d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:) d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:) d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s if (iflag_trac.eq.1) then tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:) d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s endif c print*,"d_t_ajs1=",d_t_ajs(1,:)*dtime c print*,"d_t_ajs2=",d_t_ajs(klon/2,:)*dtime c print*,"d_t_ajs3=",d_t_ajs(klon,:)*dtime c print*,"d_u_ajs=",d_u_ajs(klon/2,:)*dtime c print*,"d_v_ajs=",d_v_ajs(klon/2,:)*dtime endif c IF (if_ebil.ge.2) THEN ztit='after dry_adjust' CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(airephy,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, sens e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) END IF c==================================================================== c MICROPHYSIQUE ET CHIMIE c==================================================================== d_tr_mph(:,:,:)=0. d_tr_kim(:,:,:)=0. c on recupere tr_seri inchange, d_tr_micro, d_tr_chim, tous les trois sur nqmax c on recupere aussi qaer pour le mettre dans les sorties c si microfi=1, sortie de qaer(1:nmicro) c si nmicro != nqmax et si chimi, sortie de tr_seri(nmicro+1:nqmax) c faire un test comme pour rayonnement, avec chimi en plus comme flag, c pour voir si chimie appelee -> bouleen, qui passe dans phytrac. c faut aussi le pas de temps chimique: dtimechim, a passer.. appel_chim = 0 IF (MOD(itapchim,chimpas).EQ.0) THEN c print*,'CHIMIE ', c $ ' (itapchim=',itapchim,'/chimpas=',chimpas,')' appel_chim = 1 itapchim = 0 ENDIF itapchim = itapchim + 1 if (iflag_trac.eq.1) then c call begintime(tt0) call phytrac (debut,lafin, . nqmax,nmicro,dtime,appel_chim,dtimechim, . paprs,pplay,delp,t,rmu0,fract,pdecli,zls, . yu1,yv1,zzlev,zzlay,ftsol, . tr_seri,qaer,d_tr_mph,d_tr_kim, . fclat,reservoir) c call endtime(tt0,tt1) c ttphytra=ttphytra+tt1 c ----- ICI on ajuste radsol en tenant compte du flux de chaleur latente c d'evaporation du reservoir. c NOTE : c'est pas tres elegant mais ca permet d'eviter d'aller c toucher a clmain. if (clouds.eq.1) then radsol(:) = radsol(:)+fclat(:) !test pas de flx de chaleur latente endif if (microfi.ge.1) then tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro) . + d_tr_mph(:,:,1:nmicro)*dtime endif c PAS ELEGANT mais je n'ai pas trouve d'autres solutions : c Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs c retourne des traceurs nuls et il y a parfois des valeurs negatives qui trainent. c Pour ne diffuser le probleme, on force les valeurs negatives a ZERO. DO iq=1,nmicro DO i=1,klon DO l=1,klev if (tr_seri(i,l,iq).lt.0.) then tr_seri(i,l,iq) = 0. endif ENDDO ENDDO ENDDO c condensation: c NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!! if ((clouds.eq.1.or.(chimi)).and.nqmax.gt.nmicro) then tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax) . + d_tr_mph(:,:,nmicro+1:nqmax)*dtime endif if ((chimi).and.(nqmax.gt.nmicro)) then c chimie: tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime endif endif c ch4=0. c do l=1,llm c ch4(1,l) = tr_seri(1,l,ich4) c do j=2,jjm c ig0=1+(j-2)*iim c do i=1,iim c ch4(j,l)= ch4(j,l) + tr_seri(ig0+i,l,ich4)/iim c enddo c enddo c ch4(jjm+1,l) = tr_seri(klon,l,ich4) c enddo c do j=1,jjm+1 c write(501,*) j,ch4(j,1) c enddo c do l=1,llm c write(502,'(I3,49(ES24.17,1X))') l, c & (ch4(j,l),j=1,jjm+1) c enddo c write(501,*) "" c write(502,*) "" c------------------ c test condensation c do i=1,nqmax c if(tname(i).eq."HCN") then c print*,"HCN=" c do k=1,klev c print*,k,tr_seri(klon/2,k,i),d_tr_mph(klon/2,k,i)*dtime c v ,d_tr_kim(klon/2,k,i)*dtime c enddo c stop c endif c enddo c------------------ c==================================================================== c RAYONNEMENT c==================================================================== IF (MOD(itaprad,radpas).EQ.0) THEN c print*,'RAYONNEMENT ', c $ ' (itaprad=',itaprad,'/radpas=',radpas,')' c ATTENTION, (klon/2) ne marche pas toujours............ c print*,"radsol avant radlwsw=",radsol(klon/2) c print*,"solsw avant radlwsw=",solsw(klon/2) c print*,"sollw avant radlwsw=",sollw(klon/2) c print*,"avant radlwsw" c ---------------- c Calcul du rayon moyen des gouttes et des fractions volumique pour le TR c ---------------- IF (clouds.eq.1) THEN DO i=1,klon DO j=1,klev rmcbar(i,j)=rmcbar(i,j)/MAX(FLOAT(ncount(i,j)),1.) xfbar(i,j,:)=xfbar(i,j,:)/MAX(FLOAT(ncount(i,j)),1.) ENDDO ENDDO ENDIF c call begintime(tt0) CALL radlwsw e (dist, rmu0, fract, dtimerad, zzlev, e paprs, pplay,ftsol, t_seri, nqmax, nmicro, c tr_seri, qaer, s heat,cool,radsol, s topsw,toplw,solsw,sollw, s sollwdown, s lwnet, swnet) c call endtime(tt0,tt1) c ttrad=ttrad+tt1 c print*,"apres radlwsw" c mise a zero du rayon moyen des gouttes et des fractions volumique IF (clouds.eq.1) THEN rmcbar(:,:) = 0. xfbar(:,:,:) = 0. ncount(:,:) = 0 ENDIF c print*,"radsol apres radlwsw=",radsol(klon/2) c print*,"solsw apres radlwsw=",solsw(klon/2) c print*,"sollw apres radlwsw=",sollw(klon/2) itaprad = 0 DO k = 1, klev DO i = 1, klon dtrad(i,k) = heat(i,k)-cool(i,k) !K/s ENDDO ENDDO c print*,"heat (K/s) =",heat(klon/2,:) c print*,"cool (K/s) =",cool(klon/2,:) c print*,"dtrad1 (K/s) =",dtrad(1,:) c print*,"dtrad2 (K/s) =",dtrad(klon/2,:) c print*,"dtrad3 (K/s) =",dtrad(klon,:) ENDIF itaprad = itaprad + 1 c==================================================================== c c Ajouter la tendance des rayonnements (tous les pas) c DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime ENDDO ENDDO IF (if_ebil.ge.2) THEN ztit='after rad' CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(airephy,ztit,ip_ebil e , topsw, toplw, solsw, sollw, zero_v e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) END IF c c==================================================================== c+jld ec_conser DO k = 1, klev DO i = 1, klon d_t_ec(i,k)=0.5/cpdet(t_seri(i,k)) $ *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2) t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k) d_t_ec(i,k) = d_t_ec(i,k)/dtime END DO END DO do i=1,klon flux_ec(i,1) = 0.0 do j=2,klev flux_ec(i,j) = flux_ec(i,j-1) . +cpdet(t_seri(i,j-1))/RG*d_t_ec(i,j-1)*delp(i,j-1) enddo enddo c-jld ec_conser c==================================================================== IF (if_ebil.ge.1) THEN ztit='after physic' CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) C Comme les tendances de la physique sont ajoute dans la dynamique, C on devrait avoir que la variation d'entalpie par la dynamique C est egale a la variation de la physique au pas de temps precedent. C Donc la somme de ces 2 variations devrait etre nulle. call diagphy(airephy,ztit,ip_ebil e , topsw, toplw, solsw, sollw, sens e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) C d_h_vcol_phy=d_h_vcol C END IF C c==================================================================== c Calcul des gravity waves FLOTT c==================================================================== c c if (ok_orodr.or.ok_gw_nonoro) then cc CALCUL DE N2 c do i=1,klon c do k=2,klev c ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2. c zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1)) c enddo c enddo c call t2tpot(klon*klev,ztlev, ztetalev,zpklev) c call t2tpot(klon*klev,t_seri,ztetalay,ppk) c do i=1,klon c do k=2,klev c zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1) c zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG c zn2(i,k) = RG*zdtetalev(i,k)/(ztetalev(i,k)*zdzlev(i,k)) c zn2(i,k) = max(zn2(i,k),1.e-12) ! securite c enddo c enddo c c endif c cc ----------------------------ORODRAG c IF (ok_orodr) THEN cc cc selection des points pour lesquels le shema est actif: c igwd=0 c DO i=1,klon c itest(i)=0 cc IF ((zstd(i).gt.10.0)) THEN c IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN c itest(i)=1 c igwd=igwd+1 c idx(igwd)=i c ENDIF c ENDDO cc igwdim=MAX(1,igwd) cc cc A ADAPTER POUR TITAN !!! c CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2, c e zmea,zstd, zsig, zgam, zthe,zpic,zval, c e igwd,idx,itest, c e t_seri, u_seri, v_seri, c s zulow, zvlow, zustrdr, zvstrdr, c s d_t_oro, d_u_oro, d_v_oro) c cc ajout des tendances c t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:) c d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s c u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:) c d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s c v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:) c d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s cc c ELSE c d_t_oro = 0. c d_u_oro = 0. c d_v_oro = 0. c zustrdr = 0. c zvstrdr = 0. cc c ENDIF ! fin de test sur ok_orodr cc cc ----------------------------OROLIFT c IF (ok_orolf) THEN cc cc selection des points pour lesquels le shema est actif: c igwd=0 c DO i=1,klon c itest(i)=0 c IF ((zpic(i)-zmea(i)).GT.100.) THEN c itest(i)=1 c igwd=igwd+1 c idx(igwd)=i c ENDIF c ENDDO cc igwdim=MAX(1,igwd) cc cc A ADAPTER POUR VENUS!!! cc CALL lift_noro(klon,klev,dtime,paprs,pplay, cc e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval, cc e igwd,idx,itest, cc e t_seri, u_seri, v_seri, cc s zulow, zvlow, zustrli, zvstrli, cc s d_t_lif, d_u_lif, d_v_lif ) c cc cc ajout des tendances c t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:) c d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s c u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:) c d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s c v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:) c d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s cc c ELSE c d_t_lif = 0. c d_u_lif = 0. c d_v_lif = 0. c zustrli = 0. c zvstrli = 0. cc c ENDIF ! fin de test sur ok_orolf c cc ---------------------------- NON-ORO GRAVITY WAVES c IF(ok_gw_nonoro) then c c call flott_gwd_ran(klon,klev,dtime,pplay,zn2, c e t_seri, u_seri, v_seri, c o zustrhi,zvstrhi, c o d_t_hin, d_u_hin, d_v_hin) c cc ajout des tendances c c t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:) c d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s c u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:) c d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s c v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:) c d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s c c ELSE c d_t_hin = 0. c d_u_hin = 0. c d_v_hin = 0. c zustrhi = 0. c zvstrhi = 0. c c ENDIF ! fin de test sur ok_gw_nonoro c c==================================================================== c Transport de ballons c==================================================================== if (ballons.eq.1) then CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond, c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM) C t,pplay,u,v,zphi) ! alt above planet average radius endif !ballons c==================================================================== c Bilan de mmt angulaire c==================================================================== if (bilansmc.eq.1) then CMODDEB FLOTT C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE) C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE DO i = 1, klon zustrph(i)=0. zvstrph(i)=0. zustrcl(i)=0. zvstrcl(i)=0. ENDDO DO k = 1, klev DO i = 1, klon zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* c (paprs(i,k)-paprs(i,k+1))/rg zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* c (paprs(i,k)-paprs(i,k+1))/rg zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)* c (paprs(i,k)-paprs(i,k+1))/rg zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)* c (paprs(i,k)-paprs(i,k+1))/rg ENDDO ENDDO CALL aaam_bud (27,klon,klev,rjourvrai,gmtime, C ra,rg,romega, C rlatd,rlond,pphis, C zustrdr,zustrli,zustrcl, C zvstrdr,zvstrli,zvstrcl, C paprs,u,v) CCMODFIN FLOTT endif !bilansmc c======================================================================= c SORTIES c======================================================================= c Convertir les incrementations en tendances c DO k = 1, klev DO i = 1, klon d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime ENDDO ENDDO c print*,"vnatphy=",v(705,:) c print*,"unatphy=",u(705,:) c DO iq = 1, nqmax DO k = 1, klev DO i = 1, klon d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime ENDDO ENDDO ENDDO c------------------------ c Calcul moment cinetique c------------------------ c TEST c mangtot = 0.0 c DO k = 1, klev c DO i = 1, klon c mang(i,k) = RA*cos(rlatd(i)*RPI/180.) c . *(u_seri(i,k)+RA*cos(rlatd(i)*RPI/180.)*ROMEGA) c . *airephy(i)*delp(i,k)/RG c mangtot=mangtot+mang(i,k) c ENDDO c ENDDO c print*,"Moment cinetique total = ",mangtot c c------------------------ c c Sauvegarder les valeurs de t et u a la fin de la physique: c DO k = 1, klev DO i = 1, klon u_ancien(i,k) = u_seri(i,k) t_ancien(i,k) = t_seri(i,k) ENDDO ENDDO c c============================================================= c Ecriture des sorties c============================================================= #ifdef CPP_IOIPSL #ifdef histmth #include "write_histmth.h" #endif #ifdef histday #include "write_histday.h" #endif #ifdef histins #include "write_histins.h" #endif #endif c ---------- Sorties testphys1d ------------- if (klon.eq.1) then call writeg1d(klon,klev,t_seri,"Temp","Temperature") call writeg1d(klon,1,ftsol,"tsurf","Surface Temp") DO k = 1, klev DO i = 1, klon c tmpout(i,k) = real(heat(i,k)) tmpout(i,k) = heat(i,k) ENDDO ENDDO call writeg1d(klon,klev,tmpout, . "heat","Solar heating") DO k = 1, klev DO i = 1, klon c tmpout(i,k) = real(dtrad(i,k)) tmpout(i,k) = dtrad(i,k) ENDDO ENDDO call writeg1d(klon,klev,tmpout, . "dtrad","IR cooling") call writeg1d(klon,klev,lwnet,"lwnet","Net LW flux") call writeg1d(klon,klev,swnet,"swnet","Net SW flux") call writeg1d(klon,klev,fluxt,"flux_vdf","Turbulent flux") call writeg1d(klon,klev,flux_ajs,"flux_ajs","Dry adjust. flux") call writeg1d(klon,klev,flux_ec,"flux_ec","Ec flux") c call writeg1d(klon,1,solsw,"surfsw","Net SW flux at surface") c call writeg1d(klon,1,sollw,"surflw","Net LW flux at surface") call writeg1d(klon,1,radsol,"surfnet","Net flux at surface") call writeg1d(klon,klev,d_t_vdf,"dt_vdf","DT from clmain") call writeg1d(klon,klev,d_t_ajs,"dt_ajs","DT from ajsec") call writeg1d(klon,klev,d_t_ec,"dt_ec","DT from Ec") endif c==================================================================== c Si c'est la fin, il faut conserver l'etat de redemarrage c==================================================================== c IF (lafin) THEN itau_phy = itau_phy + itap c REMETTRE TOUS LES PARAMETRES POUR OROGW... A FAIRE POUR TITAN CALL phyredem ("restartphy.nc",dtime,radpas,chimpas, . rlatd, rlond, ftsol, ftsoil, . falbe, . solsw, sollw,dlw, . radsol,reservoir, c . zmea,zstd,zsig,zgam,zthe,zpic,zval, . t_ancien) c--------------FLOTT CMODEB LOTT C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES C DU BILAN DE MOMENT ANGULAIRE. if (bilansmc.eq.1) then write(*,*)'Fermeture de aaam_bud.out (FL Vous parle)' close(27) close(28) endif !bilansmc CMODFIN c------------- c--------------SLEBONNOIS C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES C DES BALLONS if (ballons.eq.1) then write(*,*)'Fermeture des ballons*.out' close(30) close(31) close(32) close(33) close(34) endif !ballons c------------- ENDIF RETURN END *********************************************************************** *********************************************************************** *********************************************************************** *********************************************************************** *********************************************************************** *********************************************************************** *********************************************************************** *********************************************************************** SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit) IMPLICIT none c c Tranformer une variable de la grille physique a c la grille d'ecriture c INTEGER nfield,nlon,iim,jjmp1, jjm REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield) c INTEGER i, n, ig c jjm = jjmp1 - 1 DO n = 1, nfield DO i=1,iim ecrit(i,n) = fi(1,n) ecrit(i+jjm*iim,n) = fi(nlon,n) ENDDO DO ig = 1, nlon - 2 ecrit(iim+ig,n) = fi(1+ig,n) ENDDO ENDDO RETURN END