! ! $Id: calfis_p.F 1279 2009-12-10 09:02:56Z abarral $ ! 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, $ clesphy0, $ pdufi, $ pdvfi, $ pdhfi, $ pdqfi, $ pdpsfi) #ifdef CPP_EARTH ! Ehouarn: For now, calfis_p needs Earth physics c c Auteur : P. Le Van, F. Hourdin c ......... USE dimphy USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root USE parallel, ONLY : omp_chunk, using_mpi USE mod_interface_dyn_phys USE Write_Field Use Write_field_p USE Times USE IOPHY USE infotrac 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 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" INTEGER ngridmx PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm ) #include "comconst.h" #include "comvert.h" #include "comgeom2.h" #include "control.h" #ifdef CPP_MPI include 'mpif.h' #endif c Arguments : c ----------- LOGICAL lafin REAL heure REAL pvcov(iip1,jjm,llm) REAL pucov(iip1,jjp1,llm) REAL pteta(iip1,jjp1,llm) REAL pmasse(iip1,jjp1,llm) REAL pq(iip1,jjp1,llm,nqtot) REAL pphis(iip1,jjp1) REAL pphi(iip1,jjp1,llm) c REAL pdvcov(iip1,jjm,llm) REAL pducov(iip1,jjp1,llm) REAL pdteta(iip1,jjp1,llm) REAL pdq(iip1,jjp1,llm,nqtot) c REAL pps(iip1,jjp1) REAL pp(iip1,jjp1,llmp1) REAL ppk(iip1,jjp1,llm) c REAL pdvfi(iip1,jjm,llm) REAL pdufi(iip1,jjp1,llm) REAL pdhfi(iip1,jjp1,llm) REAL pdqfi(iip1,jjp1,llm,nqtot) REAL pdpsfi(iip1,jjp1) INTEGER longcles PARAMETER ( longcles = 20 ) REAL clesphy0( longcles ) 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(:) c REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:) REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:) c 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 c REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:) REAL,ALLOCATABLE,SAVE :: zplay_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 c$OMP THREADPRIVATE(zplev_omp,zplay_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) 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 c cIM diagnostique PVteta, Amip2 INTEGER ntetaSTD PARAMETER(ntetaSTD=3) REAL rtetaSTD(ntetaSTD) DATA rtetaSTD/350., 380., 405./ REAL PVteta(klon,ntetaSTD) REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique REAL SSUM LOGICAL firstcal, debut DATA firstcal/.true./ SAVE firstcal,debut c$OMP THREADPRIVATE(firstcal,debut) REAL, intent(in):: jD_cur, jH_cur 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 c c----------------------------------------------------------------------- c c 1. Initialisations : c -------------------- c klon=klon_mpi PVteta(:,:)=0. c IF ( firstcal ) THEN debut = .TRUE. IF (ngridmx.NE.2+(jjm-1)*iim) THEN PRINT*,'STOP dans calfis' PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim' PRINT*,' ngridmx jjm iim ' PRINT*,ngridmx,jjm,iim STOP ENDIF c$OMP MASTER 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)) c$OMP END MASTER c$OMP BARRIER ELSE debut = .FALSE. ENDIF c c 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 : c 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 unskap = 1./ kappa c 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 c c c 43. temperature naturelle (en K) et pressions milieux couches . c --------------------------------------------------------------- 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 ztfi(ig0,l) = pteta(i,j,l) * pksurcp enddo ENDDO c$OMP END DO NOWAIT c 43.bis traceurs c --------------- c DO iq=1,nqtot iiq=niadv(iq) 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,iiq) enddo ENDDO c$OMP END DO NOWAIT ENDDO 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 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) THEN c cIM calcul PV a teta=350, 380, 405K CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta, $ ztfi,zplay,zplev, $ ntetaSTD,rtetaSTD,PVteta) c 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 --------------------- c$OMP BARRIER if (first_omp) then klon=klon_omp allocate(zplev_omp(klon,llm+1)) allocate(zplay_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(zdpsrf_omp(klon)) allocate(flxwfi_omp(klon,llm)) 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 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 c$OMP BARRIER if (planet_type=="earth") then #ifdef CPP_EARTH CALL physiq (klon, . llm, . debut, . lafin, . jD_cur, . jH_cur, . dtphys, . zplev_omp, . zplay_omp, . zphi_omp, . zphis_omp, . presnivs_omp, . clesphy0, . 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) #endif endif !of if (planet_type=="earth") 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