! ! $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, . flxmw, . d_u, d_v, d_t, d_qx, d_ps) c====================================================================== c c Modifications pour la physique de Venus c S. Lebonnois (LMD/CNRS) Septembre 2005 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 rjour---input-R-numero du jour de l'experience c gmtime--input-R-fraction de la journee (0 a 1) 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 flxmw---input-R-flux de masse vertical en kg/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 ! not needed; histcom is included in ioipsl USE infotrac use dimphy USE geometry_mod, only: longitude, latitude, cell_area, dx, dy USE mod_phys_lmdz_para, only : is_parallel,jj_nb USE phys_state_var_mod ! Variables sauvegardees de la physique USE write_field_phy USE iophy USE cpdet_mod, only: cpdet, t2tpot USE chemparam_mod USE conc USE compo_hedin83_mod2 use moyzon_mod, only: tmoy ! use ieee_arithmetic use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy use mod_grid_phy_lmdz, only: nbp_lon use logic_mod, only: iflag_trac IMPLICIT none c====================================================================== c CLEFS CPP POUR LES IO c ===================== c#define histhf #define histday #define histmth #define histins c====================================================================== #include "dimsoil.h" #include "clesphys.h" #include "iniprint.h" #include "timerad.h" #include "tabcontrol.h" #include "nirdata.h" #include "hedin.h" c====================================================================== LOGICAL ok_journe ! sortir le fichier journalier save ok_journe c PARAMETER (ok_journe=.true.) c LOGICAL ok_mensuel ! sortir le fichier mensuel save ok_mensuel c PARAMETER (ok_mensuel=.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 d_u_dyn(klon,klev) REAL d_t_dyn(klon,klev) REAL flxmw(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) logical ok_hf real ecrit_hf integer nid_hf save ok_hf, ecrit_hf, nid_hf #ifdef histhf data ok_hf,ecrit_hf/.true.,0.25/ #else data ok_hf/.false./ #endif c Variables propres a la physique c INTEGER,save :: itap ! compteur pour la physique REAL delp(klon,klev) ! epaisseur d'une couche REAL omega(klon,klev) ! vitesse verticale en Pa/s 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 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(:,:) 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 REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee 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 REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2) c====================================================================== c c Declaration des procedures appelees c EXTERNAL ajsec ! ajustement sec EXTERNAL clmain ! couche limite EXTERNAL hgardfou ! verifier les temperatures c 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 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 nirco2abs EXTERNAL nir_leedat EXTERNAL nltecool EXTERNAL nlte_tcool EXTERNAL nlte_setup EXTERNAL blendrad EXTERNAL nlthermeq EXTERNAL euvheat EXTERNAL param_read EXTERNAL param_read_e107 EXTERNAL conduction EXTERNAL molvis EXTERNAL moldiff_red 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 REAL tmpout(klon,klev) ! K s-1 INTEGER itaprad SAVE itaprad c REAL dist, rmu0(klon), fract(klon) REAL zdtime, zlongi c INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev c REAL zphi(klon,klev) REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2 real tsurf(klon) c va avec nlte_tcool INTEGER ierr_nlte REAL varerr c Variables du changement c c ajs: ajustement sec c vdf: couche limite (Vertical DiFfusion) 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 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 Tendencies due to radiative scheme [K/s] c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv c are not computed at each physical timestep c therefore, they are defined and saved in phys_state_var_mod c Tendencies due to molecular viscosity and conduction real d_t_conduc(klon,klev) ! [K/s] real d_u_molvis(klon,klev) ! (m/s) /s real d_v_molvis(klon,klev) ! (m/s) /s c Tendencies due to molecular diffusion real d_q_moldif(klon,klev,nqmax) 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) REAL :: d_tr(klon,klev,nqmax) c Champ de modification de la temperature par rapport a VIRAII REAL delta_temp(klon,klev) c SAVE delta_temp REAL mat_dtemp(33,50) SAVE mat_dtemp c Variables tendance sedimentation REAL :: d_tr_sed(klon,klev,2) REAL :: d_tr_ssed(klon) c c pour ioipsl INTEGER nid_day, nid_mth, nid_ins SAVE nid_day, nid_mth, nid_ins 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 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 ! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!! if (is_parallel) then bilansmc = 0 ballons = 0 endif 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(source(klon,nqmax)) 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_journe, ok_mensuel, . ok_instan, . if_ebil) call phys_state_var_init c c Initialising Hedin model for upper atm c (to be revised when coupled to chemistry) : call conc_init c c Initialiser les compteurs: c itap = 0 itaprad = 0 c c Lecture startphy.nc : c CALL phyetat0 ("startphy.nc") c dtime est defini dans tabcontrol.h et 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) CALL printflag( ok_journe,ok_instan ) c c--------- c FLOTT IF (ok_orodr) THEN DO i=1,klon rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0) ENDDO CALL SUGWD(klon,klev,paprs,pplay) DO i=1,klon zuthe(i)=0. zvthe(i)=0. if(zstd(i).gt.10.)then zuthe(i)=(1.-zgam(i))*cos(zthe(i)) zvthe(i)=(1.-zgam(i))*sin(zthe(i)) endif ENDDO 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--------- c INITIALIZE THERMOSPHERIC PARAMETERS c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (callthermos) then if(solvarmod.eq.0) call param_read if(solvarmod.eq.1) call param_read_e107 endif c Initialisation (recomputed in concentration2) do ig=1,klon do j=1,klev rnew(ig,j)=R cpnew(ig,j)=cpdet(tmoy(j)) mmean(ig,j)=RMD akknew(ig,j)=1.e-4 enddo c stop enddo IF(callthermos.or.callnlte.or.callnirco2) THEN call compo_hedin83_init2 ENDIF if (callnlte.and.nltemodel.eq.2) call nlte_setup if (callnirco2.and.nircorr.eq.1) call nir_leedat c--------- c c Verifications: c 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*REAL(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*ecriphy) ! tous les ecritphy 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*ecriphy) ! 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 histhf #include "ini_histhf.h" #endif #ifdef histday #include "ini_histday.h" #endif #ifdef histmth #include "ini_histmth.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 c--------- c Ecriture fichier initialisation c PRINT*,'Ecriture Initial_State.csv' c OPEN(88,file='Trac_Point.csv', c & form='formatted') c--------- c--------- c Initialisation des parametres des nuages c=============================================== if ((nlon .EQ. 1) .AND. ok_cloud) then PRINT*,'Open profile_cloud_parameters.csv' OPEN(66,file='profile_cloud_parameters.csv', & form='formatted') endif if ((nlon .EQ. 1) .AND. ok_sedim) then PRINT*,'Open profile_cloud_sedim.csv' OPEN(77,file='profile_cloud_sedim.csv', & form='formatted') endif if ((nlon .GT. 1) .AND. (ok_chem.OR.ok_cloud)) then c !!! DONC 3D !!! CALL chemparam_ini() endif if ((nlon .GT. 1) .AND. ok_cloud) then c !!! DONC 3D !!! CALL cloud_ini(nlon,nlev) endif c====================================================================== c Lecture du fichier DeltaT c====================================================================== c ATTENTION tout ce qui suit est pour un 48*32*50 if (ok_deltatemp) then print*,'lecture de VenusDeltaT.txt ' open(99, form = 'formatted', status = 'old', file = & 'VenusDeltaT.dat') print*,'Ouverture de VenusDeltaT.txt ' DO ilev = 1, klev read(99,'(33(1x,e13.6))') (mat_dtemp(ilat,ilev),ilat=1,33) print*,'lecture de VenusDeltaT.txt ligne:',ilev ENDDO close(99) print*,'FIN lecture de VenusDeltaT.txt ok.' DO k = 1, klev DO i = 1, klon ilat=(latitude(i)/5.625) + 17. delta_temp(i,k)=mat_dtemp(INT(ilat),k) ENDDO ENDDO endif ENDIF ! debut c====================================================================== c====================================================================== 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(cell_area,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(cell_area,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 Calcule de vitesse verticale a partir de flux de masse verticale DO k = 1, klev DO i = 1, klon omega(i,k) = RG*flxmw(i,k) / cell_area(i) END DO END DO 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 k=1,klev DO i=1,klon zzlay(i,k)=zphi(i,k)/RG ! [m] ENDDO ENDDO DO i=1,klon zzlev(i,1)=pphis(i)/RG ! [m] ENDDO DO k=2,klev DO i=1,klon z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k)) z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k)) zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2) ENDDO ENDDO DO i=1,klon zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev)) ENDDO 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 Orbite et eclairement c==================================================================== c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0. c donc pas de variations de Ls, ni de dist. c La seule chose qui compte, c'est la rotation de la planete devant c le Soleil... zlongi = 0.0 dist = 0.72 ! en UA c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite c a sa valeur, et prendre en compte leur evolution, c il faudra refaire un orbite.F... c CALL orbite(zlongi,dist) IF (cycle_diurne) THEN zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s) CALL zenang(zlongi,gmtime,zdtime,latitude,longitude,rmu0,fract) ELSE call mucorr(klon,zlongi,latitude,rmu0,fract) ENDIF c==================================================================== c Calcul des tendances traceurs c==================================================================== if (iflag_trac.eq.1) then if (tr_scheme.eq.1) then ! Case 1: pseudo-chemistry with relaxation toward fixed profile call phytrac_relax (debut,lafin,nqmax, I nlon,nlev,dtime,pplay, O tr_seri) elseif (tr_scheme.eq.2) then ! Case 2: surface emission ! For the moment, inspired from Mars version ! However, the variable 'source' could be used in physiq ! so the call to phytrac_emiss could be to initialise it. call phytrac_emiss ( (rjourvrai+gmtime)*RDAY, I debut,lafin,nqmax, I nlon,nlev,dtime,paprs, I latitude,longitude, O tr_seri) elseif (tr_scheme.eq.3) then ! identical to ok_chem.or.ok_cloud ! Case 3: Full chemistry and/or clouds if (ok_deltatemp) then ! PRINT*,'Def de delta_temp' DO k = 1, klev DO i = 1, klon ilat=(latitude(i)/5.625) + 17. ! PRINT*,INT(ilat),latitude(i),mat_dtemp(INT(ilat),k) delta_temp(i,k)=mat_dtemp(INT(ilat),k) ENDDO ENDDO endif ! Pas appel concentrations car BUG ... comprend pas ! valeur de mmean trop hautes 5,25E+29 ! call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax) if (ok_deltatemp) then ! Utilisation du champ de temperature modifie call phytrac_chimie( I debut, I gmtime, I nqmax, I klon, I latitude, I longitude, I nlev, I dtime, I t_seri+delta_temp, I pplay, O tr_seri) else call phytrac_chimie( I debut, I gmtime, I nqmax, I klon, I latitude, I longitude, I nlev, I dtime, I t_seri, I pplay, O tr_seri) endif c CALL WriteField_phy('Pression',pplay,nlev) c CALL WriteField_phy('PressionBnd',paprs,nlev+1) c CALL WriteField_phy('Temp',t_seri,nlev) c IF (ok_cloud) THEN c CALL WriteField_phy('NBRTOT',NBRTOT,nlev) c ENDIF c CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev) c CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev) if (ok_sedim) then if (ok_deltatemp) then ! Utilisation du champ de temperature modifie CALL new_cloud_sedim( I klon, I nlev, I dtime, I pplay, I paprs, I t_seri+delta_temp, I tr_seri, O d_tr_sed, O d_tr_ssed, I nqmax, O Fsedim) else CALL new_cloud_sedim( I klon, I nlev, I dtime, I pplay, I paprs, I t_seri, I tr_seri, O d_tr_sed, O d_tr_ssed, I nqmax, O Fsedim) endif DO k = 1, klev DO i = 1, klon c-------------------- c Ce test est necessaire pour eviter Xliq=NaN ! IF (ieee_is_nan(d_tr_sed(i,k,1)).OR. ! & ieee_is_nan(d_tr_sed(i,k,2))) THEN IF ((d_tr_sed(i,k,1).ne.d_tr_sed(i,k,1)).OR. & (d_tr_sed(i,k,2).ne.d_tr_sed(i,k,2))) THEN PRINT*,'sedim NaN PROBLEM' PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k) PRINT*,'lat-lon',i,'level',k,'dtime',dtime PRINT*,'F_sed',Fsedim(i,k) PRINT*,'===============================================' d_tr_sed(i,k,:)=0. ENDIF c-------------------- tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+ & d_tr_sed(i,k,1) tr_seri(i,k,i_h2oliq) = tr_seri(i,k,i_h2oliq)+ & d_tr_sed(i,k,2) d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime Fsedim(i,k) = Fsedim(i,k) / dtime ENDDO ENDDO Fsedim(:,klev+1) = 0. endif ! ok_sedim endif ! tr_scheme endif ! iflag_trac c c==================================================================== c Appeler la diffusion verticale (programme de couche limite) c==================================================================== c------------------------------- c VENUS 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 ! 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 longitude, latitude, dx, dy, 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) 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 TRACEURS if (iflag_trac.eq.1) then DO k = 1, klev DO i = 1, klon delp(i,k) = paprs(i,k)-paprs(i,k+1) ENDDO ENDDO DO iq=1, nqmax CALL cltrac(dtime,ycoefh,t_seri, s tr_seri(:,:,iq),source(:,iq), e paprs, pplay,delp, s d_tr_vdf(:,:,iq)) tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq) d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s ENDDO !nqmax endif IF (if_ebil.ge.2) THEN ztit='after clmain' CALL diagetpq(cell_area,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(cell_area,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 DO i = 1, klon ftsol(i) = ftsol(i) + d_ts(i) ENDDO 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 VENUS TEST ! tests: output tendencies ! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev) ! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev) ! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev) ! call writefield_phy('physiq_d_ts',d_ts,1) 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 . *(paprs(i,j-1)-paprs(i,j)) 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 endif ! tests: output tendencies ! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev) ! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev) ! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev) c IF (if_ebil.ge.2) THEN ztit='after dry_adjust' CALL diagetpq(cell_area,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(cell_area,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 RAYONNEMENT c==================================================================== c------------------------------------ c . Compute radiative tendencies : c------------------------------------ c==================================================================== IF (MOD(itaprad,radpas).EQ.0) THEN c==================================================================== dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s) c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas c------------------------------------ c . Compute mean mass, cp and R : c------------------------------------ if(callthermos) then call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax) endif cc!!! ADD key callhedin IF(callnlte.or.callthermos) THEN call compo_hedin83_mod(pplay,rmu0, & co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm) IF(ok_chem) then CC !! GG : Using only mayor species tracers abundances to compute NLTE heating/cooling CC Conversion [mmr] ---> [vmr] co2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co2)* & mmean(1:nlon,1:nlev)/M_tr(i_co2) covmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co)* & mmean(1:nlon,1:nlev)/M_tr(i_co) ovmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_o)* & mmean(1:nlon,1:nlev)/M_tr(i_o) n2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_n2)* & mmean(1:nlon,1:nlev)/M_tr(i_n2) ENDIF ENDIF c c NLTE cooling from CO2 emission c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF(callnlte) THEN if(nltemodel.eq.0.or.nltemodel.eq.1) then CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri, $ tr_seri, d_t_nlte) else if(nltemodel.eq.2) then CALL nlte_tcool(klon,klev,pplay*9.869e-6, $ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, $ ovmr_gcm,d_t_nlte,ierr_nlte,varerr ) if(ierr_nlte.gt.0) then write(*,*) $ 'WARNING: nlte_tcool output with error message', $ 'ierr_nlte=',ierr_nlte,'varerr=',varerr write(*,*)'I will continue anyway' endif endif ELSE d_t_nlte(:,:)=0. ENDIF c Find number of layers for LTE radiation calculations IF(callnlte .or. callnirco2) $ CALL nlthermeq(klon, klev, paprs, pplay) c c LTE radiative transfert / solar / IR matrix c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL radlwsw e (dist, rmu0, fract, zzlev, e paprs, pplay,ftsol, t_seri) c CO2 near infrared absorption c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ d_t_nirco2(:,:)=0. if (callnirco2) then call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri, . rmu0, fract, d_t_nirco2) endif c Net atmospheric radiative heating rate (K.s-1) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF(callnlte.or.callnirco2) THEN CALL blendrad(klon, klev, pplay,heat, & cool, d_t_nirco2,d_t_nlte, dtsw, dtlw) ELSE dtsw(:,:)=heat(:,:) dtlw(:,:)=-1*cool(:,:) ENDIF DO k=1,klev DO i=1,klon d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s ENDDO ENDDO cc--------------------------------------------- c EUV heating rate (K.s-1) c ~~~~~~~~~~~~~~~~~~~~~~~~ d_t_euv(:,:)=0. IF (callthermos) THEN c call euvheat(klon, klev,t_seri,paprs,pplay,zzlay, c $ rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm, c $ covmr_gcm, ovmr_gcm,d_t_euv ) call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay, $ rmu0,pdtphys,gmtime,rjourvrai, $ tr_seri, d_tr, d_t_euv ) DO k=1,klev DO ig=1,klon d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k) ENDDO ENDDO ENDIF ! callthermos c==================================================================== itaprad = 0 ENDIF ! radpas 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) + d_t_rad(i,k) * dtime ENDDO ENDDO ! CONDUCTION and MOLECULAR VISCOSITY c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ d_t_conduc(:,:)=0. d_u_molvis(:,:)=0. d_v_molvis(:,:)=0. IF (callthermos) THEN tsurf(:)=t_seri(:,1) call conduction(klon, klev,pdtphys, $ pplay,paprs,t_seri, $ tsurf,zzlev,zzlay,d_t_conduc) call molvis(klon, klev,pdtphys, $ pplay,paprs,t_seri, $ u,tsurf,zzlev,zzlay,d_u_molvis) call molvis(klon, klev, pdtphys, $ pplay,paprs,t_seri, $ v,tsurf,zzlev,zzlay,d_v_molvis) DO k=1,klev DO ig=1,klon t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K] u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s ENDDO ENDDO ENDIF ! -- MOLECULAR DIFFUSION --- d_q_moldif(:,:,:)=0 IF (callthermos .and. ok_chem) THEN call moldiff_red(klon, klev, nqmax, & pplay,paprs,t_seri, tr_seri, pdtphys, & zzlay,d_t_euv,d_t_conduc,d_q_moldif) ! --- update tendencies tracers --- DO iq = 1, nqmax DO k=1,klev DO ig=1,klon tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+ & d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]? ENDDO ENDDO ENDDO ENDIF ! callthermos & ok_chem c==================================================================== ! tests: output tendencies ! call writefield_phy('physiq_dtrad',dtrad,klev) IF (if_ebil.ge.2) THEN ztit='after rad' CALL diagetpq(cell_area,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(cell_area,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 Calcul des gravity waves FLOTT c==================================================================== c if (ok_orodr.or.ok_gw_nonoro) then c CALCUL DE N2 do i=1,klon do k=2,klev ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2. zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1)) enddo enddo call t2tpot(klon*klev,ztlev, ztetalev,zpklev) call t2tpot(klon*klev,t_seri,ztetalay,ppk) do i=1,klon do k=2,klev zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1) zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG zn2(i,k) = RG*zdtetalev(i,k)/(ztetalev(i,k)*zdzlev(i,k)) zn2(i,k) = max(zn2(i,k),1.e-12) ! securite enddo zn2(i,1) = 1.e-12 ! securite enddo endif c ----------------------------ORODRAG IF (ok_orodr) THEN c c selection des points pour lesquels le shema est actif: igwd=0 DO i=1,klon itest(i)=0 c IF ((zstd(i).gt.10.0)) THEN IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN itest(i)=1 igwd=igwd+1 idx(igwd)=i ENDIF ENDDO c igwdim=MAX(1,igwd) c c A ADAPTER POUR VENUS!!! CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2, e zmea,zstd, zsig, zgam, zthe,zpic,zval, e igwd,idx,itest, e t_seri, u_seri, v_seri, s zulow, zvlow, zustrdr, zvstrdr, s d_t_oro, d_u_oro, d_v_oro) c print*,"d_u_oro=",d_u_oro(klon/2,:) c ajout des tendances t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:) d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:) d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:) d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s c ELSE d_t_oro = 0. d_u_oro = 0. d_v_oro = 0. zustrdr = 0. zvstrdr = 0. c ENDIF ! fin de test sur ok_orodr c ! tests: output tendencies ! call writefield_phy('physiq_d_t_oro',d_t_oro,klev) ! call writefield_phy('physiq_d_u_oro',d_u_oro,klev) ! call writefield_phy('physiq_d_v_oro',d_v_oro,klev) c ----------------------------OROLIFT IF (ok_orolf) THEN print*,"ok_orolf NOT IMPLEMENTED !" stop c c selection des points pour lesquels le shema est actif: igwd=0 DO i=1,klon itest(i)=0 IF ((zpic(i)-zmea(i)).GT.100.) THEN itest(i)=1 igwd=igwd+1 idx(igwd)=i ENDIF ENDDO c igwdim=MAX(1,igwd) c c A ADAPTER POUR VENUS!!! c CALL lift_noro(klon,klev,dtime,paprs,pplay, c e latitude,zmea,zstd,zpic,zgam,zthe,zpic,zval, c e igwd,idx,itest, c e t_seri, u_seri, v_seri, c s zulow, zvlow, zustrli, zvstrli, c s d_t_lif, d_u_lif, d_v_lif ) c c ajout des tendances t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:) d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:) d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:) d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s c ELSE d_t_lif = 0. d_u_lif = 0. d_v_lif = 0. zustrli = 0. zvstrli = 0. c ENDIF ! fin de test sur ok_orolf c ---------------------------- NON-ORO GRAVITY WAVES IF(ok_gw_nonoro) then call flott_gwd_ran(klon,klev,dtime,pplay,zn2, e t_seri, u_seri, v_seri, o zustrhi,zvstrhi, o d_t_hin, d_u_hin, d_v_hin) c ajout des tendances t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:) d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:) d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:) d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s ELSE d_t_hin = 0. d_u_hin = 0. d_v_hin = 0. zustrhi = 0. zvstrhi = 0. ENDIF ! fin de test sur ok_gw_nonoro ! tests: output tendencies ! call writefield_phy('physiq_d_t_hin',d_t_hin,klev) ! call writefield_phy('physiq_d_u_hin',d_u_hin,klev) ! call writefield_phy('physiq_d_v_hin',d_v_hin,klev) c==================================================================== c Transport de ballons c==================================================================== if (ballons.eq.1) then CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,latitude,longitude, 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*RDAY, C ra,rg,romega, C latitude,longitude,pphis, C zustrdr,zustrli,zustrcl, C zvstrdr,zvstrli,zvstrcl, C paprs,u,v) CCMODFIN FLOTT endif !bilansmc c==================================================================== c==================================================================== c Calculer le transport de l'eau et de l'energie (diagnostique) c c A REVOIR POUR VENUS... c c CALL transp (paprs,ftsol, c e t_seri, q_seri, u_seri, v_seri, zphi, c s ve, vq, ue, uq) 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)*(paprs(i,j-1)-paprs(i,j)) enddo enddo c-jld ec_conser c==================================================================== IF (if_ebil.ge.1) THEN ztit='after physic' CALL diagetpq(cell_area,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(cell_area,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 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 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 VENUS... c mangtot = 0.0 c DO k = 1, klev c DO i = 1, klon c mang(i,k) = RA*cos(latitude(i)*RPI/180.) c . *(u_seri(i,k)+RA*cos(latitude(i)*RPI/180.)*ROMEGA) c . *cell_area(i)*(paprs(i,k)-paprs(i,k+1))/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 histhf #include "write_histhf.h" #endif #ifdef histday #include "write_histday.h" #endif #ifdef histmth #include "write_histmth.h" #endif #ifdef histins #include "write_histins.h" #endif #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 CALL phyredem ("restartphy.nc") 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