! ! $Id: calfis_p.F 1407 2010-07-07 10:31:52Z fairhead $ ! C C SUBROUTINE calfis_p(lafin, $ jD_cur, jH_cur, $ pucov, $ pvcov, $ pteta, $ pq, $ pmasse, $ pps, $ pp, $ ppk, $ pphis, $ pphi, $ pducov, $ pdvcov, $ pdteta, $ pdq, $ flxw, $ pdufi, $ pdvfi, $ pdhfi, $ pdqfi, $ pdpsfi) #ifdef CPP_PHYS ! Ehouarn: if using (parallelized) physics USE dimphy USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root USE mod_interface_dyn_phys ! USE IOPHY #endif USE parallel_lmdz, ONLY : omp_chunk, using_mpi, AllGather_Field USE Write_Field Use Write_field_p USE Times USE infotrac, ONLY: nqtot, niadv, tname USE control_mod, ONLY: planet_type, nsplit_phys USE cpdet_mod, only: tpot2t_p, t2tpot_p ! used only for zonal averages USE moyzon_mod IMPLICIT NONE c======================================================================= c c 1. rearrangement des tableaux et transformation c variables dynamiques > variables physiques c 2. calcul des termes physiques c 3. retransformation des tendances physiques en tendances dynamiques c c remarques: c ---------- c c - les vents sont donnes dans la physique par leurs composantes c naturelles. c - la variable thermodynamique de la physique est une variable c intensive : T c pour la dynamique on prend T * ( preff / p(l) ) **kappa c - les deux seules variables dependant de la geometrie necessaires c pour la physique sont la latitude pour le rayonnement et c l'aire de la maille quand on veut integrer une grandeur c horizontalement. c - les points de la physique sont les points scalaires de la c la dynamique; numerotation: c 1 pour le pole nord c (jjm-1)*iim pour l'interieur du domaine c ngridmx pour le pole sud c ---> ngridmx=2+(jjm-1)*iim c c Input : c ------- c ecritphy frequence d'ecriture (en jours)de histphy c pucov covariant zonal velocity c pvcov covariant meridional velocity c pteta potential temperature c pps surface pressure c pmasse masse d'air dans chaque maille c pts surface temperature (K) c callrad clef d'appel au rayonnement c c Output : c -------- c pdufi tendency for the natural zonal velocity (ms-1) c pdvfi tendency for the natural meridional velocity c pdhfi tendency for the potential temperature (K/s) c pdtsfi tendency for the surface temperature c c pdtrad radiative tendencies \ both input c pfluxrad radiative fluxes / and output c c======================================================================= c c----------------------------------------------------------------------- c c 0. Declarations : c ------------------ #include "dimensions.h" #include "paramet.h" #include "temps.h" #include "logic.h" INTEGER ngridmx PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm ) #include "comconst.h" #include "comvert.h" #include "comgeom2.h" #include "iniprint.h" #ifdef CPP_MPI include 'mpif.h' #endif c Arguments : c ----------- LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics REAL,INTENT(IN) :: jD_cur, jH_cur REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta ! commentaire SL: pdq ne sert que pour le calcul de pcvgq, ! qui lui meme ne sert a rien dans la routine telle qu'elle est ! ecrite, et que j'ai donc commente.... REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers ! NB: pdq is only used to compute pcvgq which is in fact not used... REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa) REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa) REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on dynamics grid ! tendencies (in */s) from the physics REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s) REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s) #ifdef CPP_PHYS c Local variables : c ----------------- INTEGER i,j,l,ig0,ig,iq,iiq REAL,ALLOCATABLE,SAVE :: zpsrf(:) REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:) REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:) REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:) REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:) ! ADAPTATION GCM POUR CP(T) REAL,ALLOCATABLE,SAVE :: zteta(:,:) REAL,ALLOCATABLE,SAVE :: zpk(:,:) ! Ces calculs ne servent pas. ! Si necessaire, decommenter ces variables et les calculs... ! REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:) ! REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:) c REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:) REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:) REAL,ALLOCATABLE,SAVE :: zdpsrf(:) REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:) REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:) REAL,ALLOCATABLE,SAVE :: zpk_omp(:,:) REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:) REAL,ALLOCATABLE,SAVE :: zphis_omp(:) REAL,ALLOCATABLE,SAVE :: presnivs_omp(:) REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:) REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:) REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:) REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:) REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:) REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:) REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:) REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:) REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:) REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Introduction du splitting (FH) ! Question pour Yann : ! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent ! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il ! soit allocatable (plutot par exemple que de passer une dimension ! dépendant du process en argument des routines) et que, du coup, ! le SAVE évite d'avoir à refaire l'allocation à chaque appel. ! Tu confirmes ? ! J'ai suivi le même principe pour les zdufic_omp ! Mais c'est surement bien que tu controles. ! REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:) REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:) REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:) REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:) REAL jH_cur_split,zdt_split LOGICAL debut_split,lafin_split INTEGER isplit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zpk_omp,zphi_omp,zphis_omp, c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp, c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp, c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp, c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp) LOGICAL,SAVE :: first_omp=.true. c$OMP THREADPRIVATE(first_omp) REAL zsin(iim),zcos(iim),z1(iim) REAL zsinbis(iim),zcosbis(iim),z1bis(iim) REAL unskap, pksurcp save unskap cIM diagnostique PVteta, Amip2 INTEGER,PARAMETER :: ntetaSTD=3 REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !! REAL PVteta(klon,ntetaSTD) REAL SSUM LOGICAL,SAVE :: firstcal=.true., debut=.true. c$OMP THREADPRIVATE(firstcal,debut) REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv INTEGER :: ierr #ifdef CPP_MPI INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status #else INTEGER,dimension(1,4) :: Status #endif INTEGER, dimension(4) :: Req REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:) integer :: k,kstart,kend INTEGER :: offset LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of ! For Titan only right now: ! to allow for 2D computation of microphys and chemistry LOGICAL,save :: flag_moyzon REAL,allocatable,save :: tmpvar(:,:) REAL,allocatable,save :: tmpvarp1(:,:) REAL,allocatable,save :: tmpvarbar(:) REAL,allocatable,save :: tmpvarbarp1(:) real :: zz1,zz2 c----------------------------------------------------------------------- c 1. Initialisations : c -------------------- klon=klon_mpi PVteta(:,:)=0. IF ( firstcal ) THEN debut = .TRUE. IF (ngridmx.NE.2+(jjm-1)*iim) THEN write(lunout,*) 'STOP dans calfis' write(lunout,*) & 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim' write(lunout,*) ' ngridmx jjm iim ' write(lunout,*) ngridmx,jjm,iim STOP ENDIF unskap = 1./ kappa c$OMP MASTER flag_moyzon = .false. if(moyzon_ch.or.moyzon_mu) then flag_moyzon = .true. allocate(tmpvar(iip1,llm)) allocate(tmpvarp1(iip1,llmp1)) allocate(tmpvarbar(llm)) allocate(tmpvarbarp1(llmp1)) endif ALLOCATE(zpsrf(klon)) ALLOCATE(zplev(klon,llm+1),zplay(klon,llm)) ALLOCATE(zphi(klon,llm),zphis(klon)) ALLOCATE(zufi(klon,llm), zvfi(klon,llm)) ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot)) ! ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm)) ! ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2)) ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm)) ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot)) ALLOCATE(zdpsrf(klon)) ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm)) ALLOCATE(flxwfi(klon,llm)) ! ADAPTATION GCM POUR CP(T) ALLOCATE(zteta(klon,llm)) ALLOCATE(zpk(klon,llm)) if (flag_moyzon) call moyzon_init c------------------------------------------------------------------ c moyennes globales pour les profils de pression et de temperature if(planet_type.eq."titan".or.planet_type.eq."venus") then call AllGather_Field(pp,iip1*jjp1,llmp1) call AllGather_Field(pteta,iip1*jjp1,llm) call AllGather_Field(ppk,iip1*jjp1,llm) call AllGather_Field(pphi,iip1*jjp1,llm) call AllGather_Field(pphis,iip1*jjp1,1) ALLOCATE(plevmoy(llm+1)) ALLOCATE(playmoy(llm)) ALLOCATE(tmoy(llm)) ALLOCATE(tetamoy(llm)) ALLOCATE(pkmoy(llm)) ALLOCATE(phimoy(0:llm)) ALLOCATE(zlevmoy(llm+1)) ALLOCATE(zlaymoy(llm)) plevmoy=0. do l=1,llmp1 do i=1,iip1 do j=1,jjp1 plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1) enddo enddo enddo tetamoy=0. pkmoy=0. phimoy=0. do i=1,iip1 do j=1,jjp1 phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1) enddo enddo do l=1,llm do i=1,iip1 do j=1,jjp1 tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1) pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1) phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1) enddo enddo enddo playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap call tpot2t_p(1,llm,tetamoy,tmoy,pkmoy) c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE: zlaymoy(:) = g*rad*rad/(g*rad-phimoy(:))-rad zlevmoy(1) = phimoy(0)/g DO l=2,llm zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l)) zz2=(plevmoy(l) +playmoy(l))/(plevmoy(l) -playmoy(l)) zlevmoy(l)=(zz1*zlaymoy(l-1)+zz2*zlaymoy(l))/(zz1+zz2) ENDDO zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm)) c------------------- c + lat index allocate(klat(klon)) do ig0=1,klon j=index_j(ig0) klat(ig0)=j enddo endif ! planet_type=titan c------------------------------------------------------------------ c$OMP END MASTER c$OMP BARRIER ELSE debut = .FALSE. ENDIF c----------------------------------------------------------------------- c 40. transformation des variables dynamiques en variables physiques: c --------------------------------------------------------------- c 41. pressions au sol (en Pascals) c ---------------------------------- c$OMP MASTER call start_timer(timer_physic) c$OMP END MASTER c$OMP MASTER !CDIR ON_ADB(index_i) !CDIR ON_ADB(index_j) do ig0=1,klon i=index_i(ig0) j=index_j(ig0) zpsrf(ig0)=pps(i,j) enddo c$OMP END MASTER c 42. pression intercouches et fonction d'Exner: c ----------------------------------------------------------------- c .... zplev definis aux (llm +1) interfaces des couches .... c .... zplay definis aux ( llm ) milieux des couches .... c ----------------------------------------------------------------- c ... Exner = cp * ( p(l) / preff ) ** kappa .... c print *,omp_rank,'klon--->',klon c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llmp1 !CDIR ON_ADB(index_i) !CDIR ON_ADB(index_j) do ig0=1,klon i=index_i(ig0) j=index_j(ig0) zplev( ig0,l ) = pp(i,j,l) enddo ENDDO c$OMP END DO NOWAIT if (flag_moyzon) then call AllGather_Field(pp,iip1*jjp1,llmp1) j=index_j(1) tmpvarp1(:,:) = pp(:,j,:) call moyzon(llmp1,tmpvarp1,tmpvarbarp1) zplevbar_mpi(1,:) = tmpvarbarp1 do ig0=2,klon j=index_j(ig0) if (j.ne.index_j(ig0-1)) then tmpvarp1(:,:) = pp(:,j,:) call moyzon(llmp1,tmpvarp1,tmpvarbarp1) zplevbar_mpi(ig0,:) = tmpvarbarp1 else zplevbar_mpi(ig0,:) = zplevbar_mpi(ig0-1,:) endif enddo endif ! ADAPTATION GCM POUR CP(T) c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm !CDIR ON_ADB(index_i) !CDIR ON_ADB(index_j) do ig0=1,klon i=index_i(ig0) j=index_j(ig0) zpk(ig0,l)=ppk(i,j,l) zteta(ig0,l)=pteta(i,j,l) enddo ENDDO c$OMP END DO NOWAIT if (flag_moyzon) then call AllGather_Field(pteta,iip1*jjp1,llm) call AllGather_Field(ppk,iip1*jjp1,llm) j=index_j(1) tmpvar(:,:) = pteta(:,j,:) call moyzon(llm,tmpvar,tmpvarbar) ztetabar_mpi(1,:) = tmpvarbar tmpvar(:,:) = ppk(:,j,:) call moyzon(llm,tmpvar,tmpvarbar) zpkbar_mpi(1,:) = tmpvarbar call tpot2t_p(1,llm,ztetabar_mpi(1,:),ztfibar_mpi(1,:), & zpkbar_mpi(1,:)) do ig0=2,klon j=index_j(ig0) if (j.ne.index_j(ig0-1)) then tmpvar(:,:) = pteta(:,j,:) call moyzon(llm,tmpvar,tmpvarbar) ztetabar_mpi(ig0,:) = tmpvarbar tmpvar(:,:) = ppk(:,j,:) call moyzon(llm,tmpvar,tmpvarbar) zpkbar_mpi(ig0,:) = tmpvarbar call tpot2t_p(1,llm,ztetabar_mpi(ig0,:),ztfibar_mpi(ig0,:), & zpkbar_mpi(ig0,:)) else zpkbar_mpi(ig0,:) = zpkbar_mpi(ig0-1,:) ztetabar_mpi(ig0,:) = ztetabar_mpi(ig0-1,:) ztfibar_mpi(ig0,:) = ztfibar_mpi(ig0-1,:) endif enddo endif c 43. temperature naturelle (en K) et pressions milieux couches . c --------------------------------------------------------------- ! ADAPTATION GCM POUR CP(T) call tpot2t_p(klon,llm,zteta,ztfi,zpk) c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm !CDIR ON_ADB(index_i) !CDIR ON_ADB(index_j) do ig0=1,klon i=index_i(ig0) j=index_j(ig0) pksurcp = ppk(i,j,l) / cpp zplay(ig0,l) = preff * pksurcp ** unskap enddo ENDDO c$OMP END DO NOWAIT if (flag_moyzon) then zplaybar_mpi(:,:) = preff * (zpkbar_mpi(:,:)/cpp)**unskap endif c 43.bis traceurs (tous intensifs) c --------------- DO iq=1,nqtot c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm !CDIR ON_ADB(index_i) !CDIR ON_ADB(index_j) do ig0=1,klon i=index_i(ig0) j=index_j(ig0) zqfi(ig0,l,iq) = pq(i,j,l,iq) enddo ENDDO c$OMP END DO NOWAIT ENDDO ! of DO iq=1,nqtot if (flag_moyzon) then DO iq=1,nqtot call AllGather_Field(pq(:,:,:,iq),iip1*jjp1,llm) j=index_j(1) tmpvar(:,:) = pq(:,j,:,iq) call moyzon(llm,tmpvar,tmpvarbar) zqfibar_mpi(1,:,iq) = tmpvarbar do ig0=2,klon j=index_j(ig0) if (j.ne.index_j(ig0-1)) then tmpvar(:,:) = pq(:,j,:,iq) call moyzon(llm,tmpvar,tmpvarbar) zqfibar_mpi(ig0,:,iq) = tmpvarbar else zqfibar_mpi(ig0,:,iq) = zqfibar_mpi(ig0-1,:,iq) endif enddo ENDDO ! of DO iq=1,nqtot endif c Geopotentiel calcule par rapport a la surface locale: c ----------------------------------------------------- CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi) CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis) c$OMP BARRIER c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm DO ig=1,klon zphi(ig,l)=zphi(ig,l)-zphis(ig) ENDDO ENDDO c$OMP END DO NOWAIT if (flag_moyzon) then call AllGather_Field(pphis,iip1*jjp1,1) call AllGather_Field(pphi,iip1*jjp1,llm) j=index_j(1) tmpvar(:,1) = pphis(:,j) call moyzon(1,tmpvar(:,1),tmpvarbar(1)) zphisbar_mpi(1) = tmpvarbar(1) tmpvar(:,:) = pphi(:,j,:) call moyzon(llm,tmpvar,tmpvarbar) zphibar_mpi(1,:) = tmpvarbar-zphisbar_mpi(1) do ig0=2,klon j=index_j(ig0) if (j.ne.index_j(ig0-1)) then tmpvar(:,1) = pphis(:,j) call moyzon(1,tmpvar(:,1),tmpvarbar(1)) zphisbar_mpi(ig0) = tmpvarbar(1) tmpvar(:,:) = pphi(:,j,:) call moyzon(llm,tmpvar,tmpvarbar) zphibar_mpi(ig0,:) = tmpvarbar-zphisbar_mpi(ig0) else zphisbar_mpi(ig0) = zphisbar_mpi(ig0-1) zphibar_mpi(ig0,:) = zphibar_mpi(ig0-1,:) endif enddo endif c c 45. champ u: c ------------ kstart=1 kend=klon if (is_north_pole) kstart=2 if (is_south_pole) kend=klon-1 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm !CDIR ON_ADB(index_i) !CDIR ON_ADB(index_j) !CDIR SPARSE do ig0=kstart,kend i=index_i(ig0) j=index_j(ig0) if (i==1) then zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j) $ + pucov(1,j,l)/cu(1,j) ) else zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j) $ + pucov(i,j,l)/cu(i,j) ) endif enddo ENDDO c$OMP END DO NOWAIT c 46.champ v: c ----------- c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm !CDIR ON_ADB(index_i) !CDIR ON_ADB(index_j) DO ig0=kstart,kend i=index_i(ig0) j=index_j(ig0) zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1) $ + pvcov(i,j,l)/cv(i,j) ) ENDDO ENDDO c$OMP END DO NOWAIT c 47. champs de vents aux pole nord c ------------------------------ c U = 1 / pi * integrale [ v * cos(long) * d long ] c V = 1 / pi * integrale [ v * sin(long) * d long ] if (is_north_pole) then c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1) DO i=2,iim z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1) ENDDO DO i=1,iim zcos(i) = COS(rlonv(i))*z1(i) zsin(i) = SIN(rlonv(i))*z1(i) ENDDO zufi(1,l) = SSUM(iim,zcos,1)/pi zvfi(1,l) = SSUM(iim,zsin,1)/pi ENDDO c$OMP END DO NOWAIT endif c 48. champs de vents aux pole sud: c --------------------------------- c U = 1 / pi * integrale [ v * cos(long) * d long ] c V = 1 / pi * integrale [ v * sin(long) * d long ] if (is_south_pole) then c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm) DO i=2,iim z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm) ENDDO DO i=1,iim zcos(i) = COS(rlonv(i))*z1(i) zsin(i) = SIN(rlonv(i))*z1(i) ENDDO zufi(klon,l) = SSUM(iim,zcos,1)/pi zvfi(klon,l) = SSUM(iim,zsin,1)/pi ENDDO c$OMP END DO NOWAIT endif IF (is_sequential.and.(planet_type=="earth")) THEN #ifdef CPP_EARTH ! PVtheta calls tetalevel, which is in the physics cIM calcul PV a teta=350, 380, 405K CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta, $ ztfi,zplay,zplev, $ ntetaSTD,rtetaSTD,PVteta) c #endif ENDIF c On change de grille, dynamique vers physiq, pour le flux de masse verticale CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi) c----------------------------------------------------------------------- c Appel de la physique: c --------------------- ! Appel de la physique: pose probleme quand on tourne ! SANS physique, car physiq.F est dans le repertoire phy[]... ! Il faut une cle CPP_PHYS ! Le fait que les arguments de physiq soient differents selon les planetes ! ne pose pas de probleme a priori. c$OMP BARRIER if (first_omp) then klon=klon_omp allocate(zplev_omp(klon,llm+1)) allocate(zplay_omp(klon,llm)) allocate(zpk_omp(klon,llm)) allocate(zphi_omp(klon,llm)) allocate(zphis_omp(klon)) allocate(presnivs_omp(llm)) allocate(zufi_omp(klon,llm)) allocate(zvfi_omp(klon,llm)) allocate(ztfi_omp(klon,llm)) allocate(zqfi_omp(klon,llm,nqtot)) allocate(zdufi_omp(klon,llm)) allocate(zdvfi_omp(klon,llm)) allocate(zdtfi_omp(klon,llm)) allocate(zdqfi_omp(klon,llm,nqtot)) allocate(zdufic_omp(klon,llm)) allocate(zdvfic_omp(klon,llm)) allocate(zdtfic_omp(klon,llm)) allocate(zdqfic_omp(klon,llm,nqtot)) allocate(zdpsrf_omp(klon)) allocate(flxwfi_omp(klon,llm)) if (flag_moyzon) call moyzon_init_omp(klon) first_omp=.false. endif klon=klon_omp offset=klon_omp_begin-1 do l=1,llm+1 do i=1,klon zplev_omp(i,l)=zplev(offset+i,l) enddo enddo do l=1,llm do i=1,klon zplay_omp(i,l)=zplay(offset+i,l) enddo enddo do l=1,llm do i=1,klon zpk_omp(i,l)=zpk(offset+i,l) enddo enddo do l=1,llm do i=1,klon zphi_omp(i,l)=zphi(offset+i,l) enddo enddo do i=1,klon zphis_omp(i)=zphis(offset+i) enddo do l=1,llm presnivs_omp(l)=presnivs(l) enddo do l=1,llm do i=1,klon zufi_omp(i,l)=zufi(offset+i,l) enddo enddo do l=1,llm do i=1,klon zvfi_omp(i,l)=zvfi(offset+i,l) enddo enddo do l=1,llm do i=1,klon ztfi_omp(i,l)=ztfi(offset+i,l) enddo enddo do iq=1,nqtot do l=1,llm do i=1,klon zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq) enddo enddo enddo do l=1,llm do i=1,klon zdufi_omp(i,l)=zdufi(offset+i,l) enddo enddo do l=1,llm do i=1,klon zdvfi_omp(i,l)=zdvfi(offset+i,l) enddo enddo do l=1,llm do i=1,klon zdtfi_omp(i,l)=zdtfi(offset+i,l) enddo enddo do iq=1,nqtot do l=1,llm do i=1,klon zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq) enddo enddo enddo do i=1,klon zdpsrf_omp(i)=zdpsrf(offset+i) enddo do l=1,llm do i=1,klon flxwfi_omp(i,l)=flxwfi(offset+i,l) enddo enddo if (flag_moyzon) then do l=1,llm+1 do i=1,klon zplevbar(i,l)=zplevbar_mpi(offset+i,l) enddo enddo do l=1,llm do i=1,klon zplaybar(i,l)=zplaybar_mpi(offset+i,l) enddo enddo do l=1,llm do i=1,klon zphibar(i,l)=zphibar_mpi(offset+i,l) enddo enddo do i=1,klon zphisbar(i)=zphisbar_mpi(offset+i) enddo do l=1,llm do i=1,klon ztfibar(i,l)=ztfibar_mpi(offset+i,l) enddo enddo do iq=1,nqtot do l=1,llm do i=1,klon zqfibar(i,l,iq)=zqfibar_mpi(offset+i,l,iq) enddo enddo enddo endif c$OMP BARRIER !$OMP MASTER ! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys !$OMP END MASTER zdt_split=dtphys/nsplit_phys zdufic_omp(:,:)=0. zdvfic_omp(:,:)=0. zdtfic_omp(:,:)=0. zdqfic_omp(:,:,:)=0. #ifdef CPP_PHYS do isplit=1,nsplit_phys jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys) debut_split=debut.and.isplit==1 lafin_split=lafin.and.isplit==nsplit_phys if (planet_type=="earth") then CALL physiq (klon, . llm, . debut_split, . lafin_split, . jD_cur, . jH_cur_split, . zdt_split, . zplev_omp, . zplay_omp, . zphi_omp, . zphis_omp, . presnivs_omp, . zufi_omp, . zvfi_omp, . ztfi_omp, . zqfi_omp, c#ifdef INCA . flxwfi_omp, c#endif . zdufi_omp, . zdvfi_omp, . zdtfi_omp, . zdqfi_omp, . zdpsrf_omp, cIM diagnostique PVteta, Amip2 . pducov, . PVteta) else if ( planet_type=="generic" ) then CALL physiq (klon, !! ngrid . llm, !! nlayer . nqtot, !! nq . tname, !! tracer names from dynamical core (given in infotrac) . debut_split, !! firstcall . lafin_split, !! lastcall . jD_cur, !! pday. see leapfrog_p . jH_cur_split, !! ptime "fraction of day" . zdt_split, !! ptimestep . zplev_omp, !! pplev . zplay_omp, !! pplay . zphi_omp, !! pphi . zufi_omp, !! pu . zvfi_omp, !! pv . ztfi_omp, !! pt . zqfi_omp, !! pq . flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq. . zdufi_omp, !! pdu . zdvfi_omp, !! pdv . zdtfi_omp, !! pdt . zdqfi_omp, !! pdq . zdpsrf_omp, !! pdpsrf . tracerdyn) !! tracerdyn <-- utilite ??? else if ( planet_type=="mars" ) then CALL physiq (klon, ! ngrid . llm, ! nlayer . nqtot, ! nq . debut_split, ! firstcall . lafin_split, ! lastcall . jD_cur, ! pday . jH_cur_split, ! ptime . zdt_split, ! ptimestep . zplev_omp, ! pplev . zplay_omp, ! pplay . zphi_omp, ! pphi . zufi_omp, ! pu . zvfi_omp, ! pv . ztfi_omp, ! pt . zqfi_omp, ! pq . flxwfi_omp, ! pw . zdufi_omp, ! pdu . zdvfi_omp, ! pdv . zdtfi_omp, ! pdt . zdqfi_omp, ! pdq . zdpsrf_omp, ! pdpsrf . tracerdyn) ! tracerdyn (somewhat obsolete) else if ((planet_type=="titan").or.(planet_type=="venus")) then CALL physiq (klon, . llm, . nqtot, . debut_split, . lafin_split, . jD_cur, . jH_cur_split, . zdt_split, . zplev_omp, . zplay_omp, . zpk_omp, . zphi_omp, . zphis_omp, . presnivs_omp, . zufi_omp, . zvfi_omp, . ztfi_omp, . zqfi_omp, . flxwfi_omp, . zdufi_omp, . zdvfi_omp, . zdtfi_omp, . zdqfi_omp, . zdpsrf_omp) else ! unknown "planet_type" write(lunout,*) "calfis_p: error, unknown planet_type: ", & trim(planet_type) stop endif ! planet_type zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:) zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:) zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:) zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:) enddo #endif ! of #ifdef CPP_PHYS zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys c$OMP BARRIER do l=1,llm+1 do i=1,klon zplev(offset+i,l)=zplev_omp(i,l) enddo enddo do l=1,llm do i=1,klon zplay(offset+i,l)=zplay_omp(i,l) enddo enddo do l=1,llm do i=1,klon zphi(offset+i,l)=zphi_omp(i,l) enddo enddo do i=1,klon zphis(offset+i)=zphis_omp(i) enddo do l=1,llm presnivs(l)=presnivs_omp(l) enddo do l=1,llm do i=1,klon zufi(offset+i,l)=zufi_omp(i,l) enddo enddo do l=1,llm do i=1,klon zvfi(offset+i,l)=zvfi_omp(i,l) enddo enddo do l=1,llm do i=1,klon ztfi(offset+i,l)=ztfi_omp(i,l) enddo enddo do iq=1,nqtot do l=1,llm do i=1,klon zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq) enddo enddo enddo do l=1,llm do i=1,klon zdufi(offset+i,l)=zdufi_omp(i,l) enddo enddo do l=1,llm do i=1,klon zdvfi(offset+i,l)=zdvfi_omp(i,l) enddo enddo do l=1,llm do i=1,klon zdtfi(offset+i,l)=zdtfi_omp(i,l) enddo enddo do iq=1,nqtot do l=1,llm do i=1,klon zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq) enddo enddo enddo do i=1,klon zdpsrf(offset+i)=zdpsrf_omp(i) enddo klon=klon_mpi 500 CONTINUE c$OMP BARRIER c$OMP MASTER call stop_timer(timer_physic) c$OMP END MASTER IF (using_mpi) THEN if (MPI_rank>0) then c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm du_send(1:iim,l)=zdufi(1:iim,l) dv_send(1:iim,l)=zdvfi(1:iim,l) ENDDO c$OMP END DO NOWAIT c$OMP BARRIER #ifdef CPP_MPI c$OMP MASTER !$OMP CRITICAL (MPI) call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401, & COMM_LMDZ,Req(1),ierr) call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402, & COMM_LMDZ,Req(2),ierr) !$OMP END CRITICAL (MPI) c$OMP END MASTER #endif c$OMP BARRIER endif if (MPI_rank0 .and. MPI_rank< MPI_Size-1) then call MPI_WAITALL(4,Req(1),Status,ierr) else if (MPI_rank>0) then call MPI_WAITALL(2,Req(1),Status,ierr) else if (MPI_rank