Context Navigation

calfis_p.F @ 113

Last change on this file since 113 was 108, checked in by slebonnois, 14 years ago
Sebastien Lebonnois: sponge layer et dissip horizontale.
File size: 30.1 KB

Rev	Line
[1]	1	!
	2	! $Id: calfis_p.F 1407 2010-07-07 10:31:52Z fairhead $
	3	!
	4	C
	5	C
	6	SUBROUTINE calfis_p(lafin,
	7	$ jD_cur, jH_cur,
	8	$ pucov,
	9	$ pvcov,
	10	$ pteta,
	11	$ pq,
	12	$ pmasse,
	13	$ pps,
	14	$ pp,
	15	$ ppk,
	16	$ pphis,
	17	$ pphi,
	18	$ pducov,
	19	$ pdvcov,
	20	$ pdteta,
	21	$ pdq,
	22	$ flxw,
	23	$ pdufi,
	24	$ pdvfi,
	25	$ pdhfi,
	26	$ pdqfi,
	27	$ pdpsfi)
[8]	28	#ifdef CPP_PHYS
[1]	29	! Ehouarn: For now, calfis_p needs Earth physics
	30	c
	31	c Auteur : P. Le Van, F. Hourdin
	32	c .........
	33	USE dimphy
	34	USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
	35	USE mod_interface_dyn_phys
	36	USE IOPHY
	37	#endif
	38	USE parallel, ONLY : omp_chunk, using_mpi
	39	USE Write_Field
	40	Use Write_field_p
	41	USE Times
	42	USE infotrac
	43	USE control_mod
	44
	45	IMPLICIT NONE
	46	c=======================================================================
	47	c
	48	c 1. rearrangement des tableaux et transformation
	49	c variables dynamiques > variables physiques
	50	c 2. calcul des termes physiques
	51	c 3. retransformation des tendances physiques en tendances dynamiques
	52	c
	53	c remarques:
	54	c ----------
	55	c
	56	c - les vents sont donnes dans la physique par leurs composantes
	57	c naturelles.
	58	c - la variable thermodynamique de la physique est une variable
	59	c intensive : T
	60	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
	61	c - les deux seules variables dependant de la geometrie necessaires
	62	c pour la physique sont la latitude pour le rayonnement et
	63	c l'aire de la maille quand on veut integrer une grandeur
	64	c horizontalement.
	65	c - les points de la physique sont les points scalaires de la
	66	c la dynamique; numerotation:
	67	c 1 pour le pole nord
	68	c (jjm-1)*iim pour l'interieur du domaine
	69	c ngridmx pour le pole sud
	70	c ---> ngridmx=2+(jjm-1)*iim
	71	c
	72	c Input :
	73	c -------
	74	c ecritphy frequence d'ecriture (en jours)de histphy
	75	c pucov covariant zonal velocity
	76	c pvcov covariant meridional velocity
	77	c pteta potential temperature
	78	c pps surface pressure
	79	c pmasse masse d'air dans chaque maille
	80	c pts surface temperature (K)
	81	c callrad clef d'appel au rayonnement
	82	c
	83	c Output :
	84	c --------
	85	c pdufi tendency for the natural zonal velocity (ms-1)
	86	c pdvfi tendency for the natural meridional velocity
[8]	87	c pdhfi tendency for the potential temperature (K/s)
[1]	88	c pdtsfi tendency for the surface temperature
	89	c
	90	c pdtrad radiative tendencies \ both input
	91	c pfluxrad radiative fluxes / and output
	92	c
	93	c=======================================================================
	94	c
	95	c-----------------------------------------------------------------------
	96	c
	97	c 0. Declarations :
	98	c ------------------
	99
	100	#include "dimensions.h"
	101	#include "paramet.h"
	102	#include "temps.h"
	103
	104	INTEGER ngridmx
	105	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
	106
	107	#include "comconst.h"
	108	#include "comvert.h"
	109	#include "comgeom2.h"
	110	#include "iniprint.h"
	111	#ifdef CPP_MPI
	112	include 'mpif.h'
	113	#endif
	114	c Arguments :
	115	c -----------
	116	LOGICAL lafin
	117	! REAL heure
	118	REAL, intent(in):: jD_cur, jH_cur
	119	REAL pvcov(iip1,jjm,llm)
	120	REAL pucov(iip1,jjp1,llm)
	121	REAL pteta(iip1,jjp1,llm)
	122	REAL pmasse(iip1,jjp1,llm)
	123	REAL pq(iip1,jjp1,llm,nqtot)
	124	REAL pphis(iip1,jjp1)
	125	REAL pphi(iip1,jjp1,llm)
	126	c
	127	REAL pdvcov(iip1,jjm,llm)
	128	REAL pducov(iip1,jjp1,llm)
	129	REAL pdteta(iip1,jjp1,llm)
[8]	130	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
	131	! qui lui meme ne sert a rien dans la routine telle qu'elle est
	132	! ecrite, et que j'ai donc commente....
[1]	133	REAL pdq(iip1,jjp1,llm,nqtot)
	134	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
	135	c
	136	REAL pps(iip1,jjp1)
	137	REAL pp(iip1,jjp1,llmp1)
	138	REAL ppk(iip1,jjp1,llm)
	139	c
	140	REAL pdvfi(iip1,jjm,llm)
	141	REAL pdufi(iip1,jjp1,llm)
	142	REAL pdhfi(iip1,jjp1,llm)
	143	REAL pdqfi(iip1,jjp1,llm,nqtot)
	144	REAL pdpsfi(iip1,jjp1)
	145
[8]	146	#ifdef CPP_PHYS
[1]	147	c Local variables :
	148	c -----------------
	149
	150	INTEGER i,j,l,ig0,ig,iq,iiq
	151	REAL,ALLOCATABLE,SAVE :: zpsrf(:)
	152	REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
	153	REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
	154	c
	155	REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
	156	REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
[8]	157	! ADAPTATION GCM POUR CP(T)
	158	REAL,ALLOCATABLE,SAVE :: zteta(:,:)
	159	REAL,ALLOCATABLE,SAVE :: zpk(:,:)
[1]	160	c
[8]	161	! Ces calculs ne servent pas.
	162	! Si necessaire, decommenter ces variables et les calculs...
	163	! REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
	164	! REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
[1]	165	c
	166	REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
	167	REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
	168	REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
	169	REAL,SAVE,ALLOCATABLE :: flxwfi(:,:) ! Flux de masse verticale sur la grille physiq
	170
	171	c
	172	REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
	173	REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
	174	REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
	175	REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
	176	REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
	177	REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:)
	178	REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
	179	REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
	180	REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
	181	REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
	182	REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
	183	REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
	184	REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
	185	REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
	186	REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq
	187
	188	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	189	! Introduction du splitting (FH)
	190	! Question pour Yann :
	191	! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
	192	! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
	193	! soit allocatable (plutot par exemple que de passer une dimension
	194	! dépendant du process en argument des routines) et que, du coup,
	195	! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
	196	! Tu confirmes ?
	197	! J'ai suivi le même principe pour les zdufic_omp
	198	! Mais c'est surement bien que tu controles.
	199	!
	200
	201	REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
	202	REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
	203	REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
	204	REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
	205	REAL jH_cur_split,zdt_split
	206	LOGICAL debut_split,lafin_split
	207	INTEGER isplit
	208	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	209
	210	c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
	211	c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
	212	c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp,
	213	c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
	214	c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)
	215
	216	LOGICAL,SAVE :: first_omp=.true.
	217	c$OMP THREADPRIVATE(first_omp)
	218
	219	REAL zsin(iim),zcos(iim),z1(iim)
	220	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
	221	REAL unskap, pksurcp
	222	c
	223	cIM diagnostique PVteta, Amip2
	224	INTEGER ntetaSTD
	225	PARAMETER(ntetaSTD=3)
	226	REAL rtetaSTD(ntetaSTD)
	227	DATA rtetaSTD/350., 380., 405./
	228	REAL PVteta(klon,ntetaSTD)
	229
	230
	231	REAL SSUM
	232
	233	LOGICAL firstcal, debut
	234	DATA firstcal/.true./
	235	SAVE firstcal,debut
	236	c$OMP THREADPRIVATE(firstcal,debut)
	237
	238	REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
	239	INTEGER :: ierr
	240	#ifdef CPP_MPI
	241	INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
	242	#else
	243	INTEGER,dimension(1,4) :: Status
	244	#endif
	245	INTEGER, dimension(4) :: Req
	246	REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
	247	integer :: k,kstart,kend
	248	INTEGER :: offset
	249	c
	250	c-----------------------------------------------------------------------
	251	c
	252	c 1. Initialisations :
	253	c --------------------
	254	c
	255
	256	klon=klon_mpi
	257
	258	PVteta(:,:)=0.
	259
	260	c
	261	IF ( firstcal ) THEN
	262	debut = .TRUE.
	263	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
	264	write(lunout,*) 'STOP dans calfis'
	265	write(lunout,*)
	266	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
	267	write(lunout,*) ' ngridmx jjm iim '
	268	write(lunout,*) ngridmx,jjm,iim
	269	STOP
	270	ENDIF
	271	c$OMP MASTER
	272	ALLOCATE(zpsrf(klon))
	273	ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
	274	ALLOCATE(zphi(klon,llm),zphis(klon))
	275	ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
	276	ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
[8]	277	! ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
	278	! ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
[1]	279	ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
	280	ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
	281	ALLOCATE(zdpsrf(klon))
	282	ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
	283	ALLOCATE(flxwfi(klon,llm))
[8]	284	! ADAPTATION GCM POUR CP(T)
[52]	285	ALLOCATE(zteta(klon,llm))
[8]	286	ALLOCATE(zpk(klon,llm))
[1]	287	c$OMP END MASTER
	288	c$OMP BARRIER
	289	ELSE
	290	debut = .FALSE.
	291	ENDIF
	292
	293	c
	294	c
	295	c-----------------------------------------------------------------------
	296	c 40. transformation des variables dynamiques en variables physiques:
	297	c ---------------------------------------------------------------
	298
	299	c 41. pressions au sol (en Pascals)
	300	c ----------------------------------
	301
	302	c$OMP MASTER
	303	call start_timer(timer_physic)
	304	c$OMP END MASTER
	305
	306	c$OMP MASTER
	307	!CDIR ON_ADB(index_i)
	308	!CDIR ON_ADB(index_j)
	309	do ig0=1,klon
	310	i=index_i(ig0)
	311	j=index_j(ig0)
	312	zpsrf(ig0)=pps(i,j)
	313	enddo
	314	c$OMP END MASTER
	315
	316
[8]	317	c 42. pression intercouches et fonction d'Exner:
[1]	318	c
	319	c -----------------------------------------------------------------
	320	c .... zplev definis aux (llm +1) interfaces des couches ....
	321	c .... zplay definis aux ( llm ) milieux des couches ....
	322	c -----------------------------------------------------------------
	323
	324	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
	325	c
	326	unskap = 1./ kappa
	327	c
	328	c print *,omp_rank,'klon--->',klon
	329	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	330	DO l = 1, llmp1
	331	!CDIR ON_ADB(index_i)
	332	!CDIR ON_ADB(index_j)
	333	do ig0=1,klon
	334	i=index_i(ig0)
	335	j=index_j(ig0)
	336	zplev( ig0,l ) = pp(i,j,l)
	337	enddo
	338	ENDDO
	339	c$OMP END DO NOWAIT
[8]	340	! ADAPTATION GCM POUR CP(T)
	341	DO l=1,llm
	342	!CDIR ON_ADB(index_i)
	343	!CDIR ON_ADB(index_j)
	344	do ig0=1,klon
	345	i=index_i(ig0)
	346	j=index_j(ig0)
	347	zpk(ig0,l)=ppk(i,j,l)
	348	zteta(ig0,l)=pteta(i,j,l)
	349	enddo
	350	ENDDO
	351	c$OMP END DO NOWAIT
	352
[1]	353	c
	354	c
	355
	356	c 43. temperature naturelle (en K) et pressions milieux couches .
	357	c ---------------------------------------------------------------
[8]	358
	359	! ADAPTATION GCM POUR CP(T)
	360	call tpot2t_p(ngridmx*llm,zteta,ztfi,zpk)
	361
[1]	362	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	363	DO l=1,llm
	364	!CDIR ON_ADB(index_i)
	365	!CDIR ON_ADB(index_j)
	366	do ig0=1,klon
	367	i=index_i(ig0)
	368	j=index_j(ig0)
	369	pksurcp = ppk(i,j,l) / cpp
	370	zplay(ig0,l) = preff * pksurcp ** unskap
[8]	371	! ztfi(ig0,l) = pteta(i,j,l) * pksurcp
[1]	372	enddo
	373
	374	ENDDO
	375	c$OMP END DO NOWAIT
	376
[8]	377	c 43.bis traceurs (tous intensifs)
[1]	378	c ---------------
	379	c
	380
	381	DO iq=1,nqtot
	382	iiq=niadv(iq)
	383	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	384	DO l=1,llm
	385	!CDIR ON_ADB(index_i)
	386	!CDIR ON_ADB(index_j)
	387	do ig0=1,klon
	388	i=index_i(ig0)
	389	j=index_j(ig0)
	390	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
	391	enddo
	392	ENDDO
	393	c$OMP END DO NOWAIT
[8]	394	ENDDO ! of DO iq=1,nqtot
[1]	395
	396
	397	c Geopotentiel calcule par rapport a la surface locale:
	398	c -----------------------------------------------------
	399
	400	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
	401
	402	CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
	403
	404	c$OMP BARRIER
	405
	406	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	407	DO l=1,llm
	408	DO ig=1,klon
	409	zphi(ig,l)=zphi(ig,l)-zphis(ig)
	410	ENDDO
	411	ENDDO
	412	c$OMP END DO NOWAIT
	413
	414
	415	c
	416	c 45. champ u:
	417	c ------------
	418
	419	kstart=1
	420	kend=klon
	421
	422	if (is_north_pole) kstart=2
	423	if (is_south_pole) kend=klon-1
	424
	425	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	426	DO l=1,llm
	427	!CDIR ON_ADB(index_i)
	428	!CDIR ON_ADB(index_j)
	429	!CDIR SPARSE
	430	do ig0=kstart,kend
	431	i=index_i(ig0)
	432	j=index_j(ig0)
	433	if (i==1) then
	434	zufi(ig0,l)= 0.5 *( pucov(iim,j,l)/cu(iim,j)
	435	$ + pucov(1,j,l)/cu(1,j) )
	436	else
	437	zufi(ig0,l)= 0.5*( pucov(i-1,j,l)/cu(i-1,j)
	438	$ + pucov(i,j,l)/cu(i,j) )
	439	endif
	440	enddo
	441	ENDDO
	442	c$OMP END DO NOWAIT
	443
	444	c 46.champ v:
	445	c -----------
	446
	447	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	448	DO l=1,llm
	449	!CDIR ON_ADB(index_i)
	450	!CDIR ON_ADB(index_j)
	451	DO ig0=kstart,kend
	452	i=index_i(ig0)
	453	j=index_j(ig0)
	454	zvfi(ig0,l)= 0.5 *( pvcov(i,j-1,l)/cv(i,j-1)
	455	$ + pvcov(i,j,l)/cv(i,j) )
	456
	457	ENDDO
	458	ENDDO
	459	c$OMP END DO NOWAIT
	460
	461	c 47. champs de vents aux pole nord
	462	c ------------------------------
	463	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	464	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	465
	466	if (is_north_pole) then
	467	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	468	DO l=1,llm
	469
	470	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
	471	DO i=2,iim
	472	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
	473	ENDDO
	474
	475	DO i=1,iim
	476	zcos(i) = COS(rlonv(i))*z1(i)
	477	zsin(i) = SIN(rlonv(i))*z1(i)
	478	ENDDO
	479
	480	zufi(1,l) = SSUM(iim,zcos,1)/pi
	481	zvfi(1,l) = SSUM(iim,zsin,1)/pi
	482
	483	ENDDO
	484	c$OMP END DO NOWAIT
	485	endif
	486
	487
	488	c 48. champs de vents aux pole sud:
	489	c ---------------------------------
	490	c U = 1 / pi * integrale [ v * cos(long) * d long ]
	491	c V = 1 / pi * integrale [ v * sin(long) * d long ]
	492
	493	if (is_south_pole) then
	494	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	495	DO l=1,llm
	496
	497	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
	498	DO i=2,iim
	499	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
	500	ENDDO
	501
	502	DO i=1,iim
	503	zcos(i) = COS(rlonv(i))*z1(i)
	504	zsin(i) = SIN(rlonv(i))*z1(i)
	505	ENDDO
	506
	507	zufi(klon,l) = SSUM(iim,zcos,1)/pi
	508	zvfi(klon,l) = SSUM(iim,zsin,1)/pi
	509	ENDDO
	510	c$OMP END DO NOWAIT
	511	endif
	512
	513
	514	IF (is_sequential) THEN
	515	c
	516	cIM calcul PV a teta=350, 380, 405K
	517	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
	518	$ ztfi,zplay,zplev,
	519	$ ntetaSTD,rtetaSTD,PVteta)
	520	c
	521	ENDIF
	522
	523	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
	524	CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
	525
	526	c-----------------------------------------------------------------------
	527	c Appel de la physique:
	528	c ---------------------
	529
[8]	530	! Appel de la physique: pose probleme quand on tourne
	531	! SANS physique, car physiq.F est dans le repertoire phy[]...
	532	! Il faut une cle CPP_PHYS
[1]	533
[8]	534	! Le fait que les arguments de physiq soient differents selon les planetes
	535	! ne pose pas de probleme a priori.
	536
	537
[1]	538	c$OMP BARRIER
	539	if (first_omp) then
	540	klon=klon_omp
	541
	542	allocate(zplev_omp(klon,llm+1))
	543	allocate(zplay_omp(klon,llm))
	544	allocate(zphi_omp(klon,llm))
	545	allocate(zphis_omp(klon))
	546	allocate(presnivs_omp(llm))
	547	allocate(zufi_omp(klon,llm))
	548	allocate(zvfi_omp(klon,llm))
	549	allocate(ztfi_omp(klon,llm))
	550	allocate(zqfi_omp(klon,llm,nqtot))
	551	allocate(zdufi_omp(klon,llm))
	552	allocate(zdvfi_omp(klon,llm))
	553	allocate(zdtfi_omp(klon,llm))
	554	allocate(zdqfi_omp(klon,llm,nqtot))
	555	allocate(zdufic_omp(klon,llm))
	556	allocate(zdvfic_omp(klon,llm))
	557	allocate(zdtfic_omp(klon,llm))
	558	allocate(zdqfic_omp(klon,llm,nqtot))
	559	allocate(zdpsrf_omp(klon))
	560	allocate(flxwfi_omp(klon,llm))
	561	first_omp=.false.
	562	endif
	563
	564
	565	klon=klon_omp
	566	offset=klon_omp_begin-1
	567
	568	do l=1,llm+1
	569	do i=1,klon
	570	zplev_omp(i,l)=zplev(offset+i,l)
	571	enddo
	572	enddo
	573
	574	do l=1,llm
	575	do i=1,klon
	576	zplay_omp(i,l)=zplay(offset+i,l)
	577	enddo
	578	enddo
	579
	580	do l=1,llm
	581	do i=1,klon
	582	zphi_omp(i,l)=zphi(offset+i,l)
	583	enddo
	584	enddo
	585
	586	do i=1,klon
	587	zphis_omp(i)=zphis(offset+i)
	588	enddo
	589
	590
	591	do l=1,llm
	592	presnivs_omp(l)=presnivs(l)
	593	enddo
	594
	595	do l=1,llm
	596	do i=1,klon
	597	zufi_omp(i,l)=zufi(offset+i,l)
	598	enddo
	599	enddo
	600
	601	do l=1,llm
	602	do i=1,klon
	603	zvfi_omp(i,l)=zvfi(offset+i,l)
	604	enddo
	605	enddo
	606
	607	do l=1,llm
	608	do i=1,klon
	609	ztfi_omp(i,l)=ztfi(offset+i,l)
	610	enddo
	611	enddo
	612
	613	do iq=1,nqtot
	614	do l=1,llm
	615	do i=1,klon
	616	zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
	617	enddo
	618	enddo
	619	enddo
	620
	621	do l=1,llm
	622	do i=1,klon
	623	zdufi_omp(i,l)=zdufi(offset+i,l)
	624	enddo
	625	enddo
	626
	627	do l=1,llm
	628	do i=1,klon
	629	zdvfi_omp(i,l)=zdvfi(offset+i,l)
	630	enddo
	631	enddo
	632
	633	do l=1,llm
	634	do i=1,klon
	635	zdtfi_omp(i,l)=zdtfi(offset+i,l)
	636	enddo
	637	enddo
	638
	639	do iq=1,nqtot
	640	do l=1,llm
	641	do i=1,klon
	642	zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
	643	enddo
	644	enddo
	645	enddo
	646
	647	do i=1,klon
	648	zdpsrf_omp(i)=zdpsrf(offset+i)
	649	enddo
	650
	651	do l=1,llm
	652	do i=1,klon
	653	flxwfi_omp(i,l)=flxwfi(offset+i,l)
	654	enddo
	655	enddo
	656
	657	c$OMP BARRIER
	658
	659	!$OMP MASTER
	660	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
	661	!$OMP END MASTER
	662	zdt_split=dtphys/nsplit_phys
	663	zdufic_omp(:,:)=0.
	664	zdvfic_omp(:,:)=0.
	665	zdtfic_omp(:,:)=0.
	666	zdqfic_omp(:,:,:)=0.
	667
	668	do isplit=1,nsplit_phys
	669
	670	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
	671	debut_split=debut.and.isplit==1
	672	lafin_split=lafin.and.isplit==nsplit_phys
	673
	674
[8]	675	if (planet_type=="earth") then
	676	CALL physiq (klon,
[1]	677	. llm,
	678	. debut_split,
	679	. lafin_split,
	680	. jD_cur,
	681	. jH_cur_split,
	682	. zdt_split,
	683	. zplev_omp,
	684	. zplay_omp,
	685	. zphi_omp,
	686	. zphis_omp,
	687	. presnivs_omp,
	688	. zufi_omp,
	689	. zvfi_omp,
	690	. ztfi_omp,
	691	. zqfi_omp,
	692	c#ifdef INCA
	693	. flxwfi_omp,
	694	c#endif
	695	. zdufi_omp,
	696	. zdvfi_omp,
	697	. zdtfi_omp,
	698	. zdqfi_omp,
	699	. zdpsrf_omp,
	700	cIM diagnostique PVteta, Amip2
	701	. pducov,
	702	. PVteta)
[8]	703	else ! a moduler pour Mars
	704	CALL physiq (klon,
	705	. llm,
	706	. debut_split,
	707	. lafin_split,
	708	. jD_cur,
	709	. jH_cur_split,
	710	. zdt_split,
	711	. zplev_omp,
	712	. zplay_omp,
	713	. zphi_omp,
	714	. zphis_omp,
	715	. presnivs_omp,
	716	. zufi_omp,
	717	. zvfi_omp,
	718	. ztfi_omp,
	719	. zqfi_omp,
	720	c#ifdef INCA
	721	. flxwfi_omp,
	722	c#endif
	723	. zdufi_omp,
	724	. zdvfi_omp,
	725	. zdtfi_omp,
	726	. zdqfi_omp,
	727	. zdpsrf_omp,
	728	cIM diagnostique PVteta, Amip2
	729	. pducov,
	730	. PVteta)
	731	endif ! planet_type
[1]	732	zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
	733	zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
	734	ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
	735	zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
	736
	737	zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
	738	zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
	739	zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
	740	zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
	741
	742	enddo
	743
	744	zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
	745	zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
	746	zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
	747	zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
	748
	749	c$OMP BARRIER
	750
	751	do l=1,llm+1
	752	do i=1,klon
	753	zplev(offset+i,l)=zplev_omp(i,l)
	754	enddo
	755	enddo
	756
	757	do l=1,llm
	758	do i=1,klon
	759	zplay(offset+i,l)=zplay_omp(i,l)
	760	enddo
	761	enddo
	762
	763	do l=1,llm
	764	do i=1,klon
	765	zphi(offset+i,l)=zphi_omp(i,l)
	766	enddo
	767	enddo
	768
	769
	770	do i=1,klon
	771	zphis(offset+i)=zphis_omp(i)
	772	enddo
	773
	774
	775	do l=1,llm
	776	presnivs(l)=presnivs_omp(l)
	777	enddo
	778
	779	do l=1,llm
	780	do i=1,klon
	781	zufi(offset+i,l)=zufi_omp(i,l)
	782	enddo
	783	enddo
	784
	785	do l=1,llm
	786	do i=1,klon
	787	zvfi(offset+i,l)=zvfi_omp(i,l)
	788	enddo
	789	enddo
	790
	791	do l=1,llm
	792	do i=1,klon
	793	ztfi(offset+i,l)=ztfi_omp(i,l)
	794	enddo
	795	enddo
	796
	797	do iq=1,nqtot
	798	do l=1,llm
	799	do i=1,klon
	800	zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
	801	enddo
	802	enddo
	803	enddo
	804
	805	do l=1,llm
	806	do i=1,klon
	807	zdufi(offset+i,l)=zdufi_omp(i,l)
	808	enddo
	809	enddo
	810
	811	do l=1,llm
	812	do i=1,klon
	813	zdvfi(offset+i,l)=zdvfi_omp(i,l)
	814	enddo
	815	enddo
	816
	817	do l=1,llm
	818	do i=1,klon
	819	zdtfi(offset+i,l)=zdtfi_omp(i,l)
	820	enddo
	821	enddo
	822
	823	do iq=1,nqtot
	824	do l=1,llm
	825	do i=1,klon
	826	zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
	827	enddo
	828	enddo
	829	enddo
	830
	831	do i=1,klon
	832	zdpsrf(offset+i)=zdpsrf_omp(i)
	833	enddo
	834
	835
	836	klon=klon_mpi
	837	500 CONTINUE
	838	c$OMP BARRIER
	839
	840	c$OMP MASTER
	841	call stop_timer(timer_physic)
	842	c$OMP END MASTER
	843
	844	IF (using_mpi) THEN
	845
	846	if (MPI_rank>0) then
	847
	848	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	849	DO l=1,llm
	850	du_send(1:iim,l)=zdufi(1:iim,l)
	851	dv_send(1:iim,l)=zdvfi(1:iim,l)
	852	ENDDO
	853	c$OMP END DO NOWAIT
	854
	855	c$OMP BARRIER
	856	#ifdef CPP_MPI
	857	c$OMP MASTER
	858	!$OMP CRITICAL (MPI)
	859	call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
	860	& COMM_LMDZ,Req(1),ierr)
	861	call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
	862	& COMM_LMDZ,Req(2),ierr)
	863	!$OMP END CRITICAL (MPI)
	864	c$OMP END MASTER
	865	#endif
	866	c$OMP BARRIER
	867
	868	endif
	869
	870	if (MPI_rank<MPI_Size-1) then
	871	c$OMP BARRIER
	872	#ifdef CPP_MPI
	873	c$OMP MASTER
	874	!$OMP CRITICAL (MPI)
	875	call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
	876	& COMM_LMDZ,Req(3),ierr)
	877	call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
	878	& COMM_LMDZ,Req(4),ierr)
	879	!$OMP END CRITICAL (MPI)
	880	c$OMP END MASTER
	881	#endif
	882	endif
	883
	884	c$OMP BARRIER
	885
	886
	887	#ifdef CPP_MPI
	888	c$OMP MASTER
	889	!$OMP CRITICAL (MPI)
	890	if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
	891	call MPI_WAITALL(4,Req(1),Status,ierr)
	892	else if (MPI_rank>0) then
	893	call MPI_WAITALL(2,Req(1),Status,ierr)
	894	else if (MPI_rank <MPI_Size-1) then
	895	call MPI_WAITALL(2,Req(3),Status,ierr)
	896	endif
	897	!$OMP END CRITICAL (MPI)
	898	c$OMP END MASTER
	899	#endif
	900
	901	c$OMP BARRIER
	902
	903	ENDIF ! using_mpi
	904
	905
	906	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	907	DO l=1,llm
	908
	909	zdufi2(1:klon,l)=zdufi(1:klon,l)
	910	zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
	911
	912	zdvfi2(1:klon,l)=zdvfi(1:klon,l)
	913	zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l)
	914
	915	pdhfi(:,jj_begin,l)=0
	916	pdqfi(:,jj_begin,l,:)=0
	917	pdufi(:,jj_begin,l)=0
	918	pdvfi(:,jj_begin,l)=0
	919
	920	if (.not. is_south_pole) then
	921	pdhfi(:,jj_end,l)=0
	922	pdqfi(:,jj_end,l,:)=0
	923	pdufi(:,jj_end,l)=0
	924	pdvfi(:,jj_end,l)=0
	925	endif
	926
	927	ENDDO
	928	c$OMP END DO NOWAIT
	929
	930	c$OMP MASTER
	931	pdpsfi(:,jj_begin)=0
	932	if (.not. is_south_pole) then
	933	pdpsfi(:,jj_end)=0
	934	endif
	935	c$OMP END MASTER
	936	c-----------------------------------------------------------------------
	937	c transformation des tendances physiques en tendances dynamiques:
	938	c ---------------------------------------------------------------
	939
	940	c tendance sur la pression :
	941	c -----------------------------------
	942	CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
	943	c
	944	c 62. enthalpie potentielle
	945	c ---------------------
[8]	946
[1]	947	kstart=1
	948	kend=klon
	949
	950	if (is_north_pole) kstart=2
	951	if (is_south_pole) kend=klon-1
[8]	952
	953	! ADAPTATION GCM POUR CP(T)
	954	call t2tpot_p(ngridmx,llm,ztfi,zteta,zpk)
[1]	955
[8]	956
[1]	957	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	958	DO l=1,llm
	959	!CDIR ON_ADB(index_i)
	960	!CDIR ON_ADB(index_j)
	961	!cdir NODEP
	962	do ig0=kstart,kend
	963	i=index_i(ig0)
	964	j=index_j(ig0)
[8]	965	! pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
	966	pdhfi(i,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
	967	! if (i==1) pdhfi(iip1,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
	968	if (i==1) then
	969	pdhfi(iip1,j,l) = (zteta(ig0,l) - pteta(i,j,l))/dtphys
	970	endif
[1]	971	enddo
	972
	973	if (is_north_pole) then
	974	DO i=1,iip1
[8]	975	! pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
	976	pdhfi(i,1,l) = (zteta(1,l) - pteta(i,1,l))/dtphys
[1]	977	enddo
	978	endif
	979
	980	if (is_south_pole) then
	981	DO i=1,iip1
[8]	982	! pdhfi(i,jjp1,l) = cpp * zdtfi(klon,l)/ ppk(i,jjp1,l)
	983	pdhfi(i,jjp1,l) = (zteta(klon,l) - pteta(i,jjp1,l))/dtphys
[1]	984	ENDDO
	985	endif
	986	ENDDO
	987	c$OMP END DO NOWAIT
	988
	989	c 62. humidite specifique
	990	c ---------------------
	991	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
	992	! DO iq=1,nqtot
	993	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	994	! DO l=1,llm
	995	!!!cdir NODEP
	996	! do ig0=kstart,kend
	997	! i=index_i(ig0)
	998	! j=index_j(ig0)
	999	! pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq)
	1000	! if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq)
	1001	! enddo
	1002	!
	1003	! if (is_north_pole) then
	1004	! do i=1,iip1
	1005	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
	1006	! enddo
	1007	! endif
	1008	!
	1009	! if (is_south_pole) then
	1010	! do i=1,iip1
	1011	! pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq)
	1012	! enddo
	1013	! endif
	1014	! ENDDO
	1015	!c$OMP END DO NOWAIT
	1016	! ENDDO
	1017
[8]	1018	c 63. traceurs (tous en intensifs)
[1]	1019	c ------------
	1020	C initialisation des tendances
	1021
	1022	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1023	DO l=1,llm
	1024	pdqfi(:,:,l,:)=0.
	1025	ENDDO
	1026	c$OMP END DO NOWAIT
	1027
	1028	C
	1029	!cdir NODEP
	1030	DO iq=1,nqtot
	1031	iiq=niadv(iq)
	1032	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1033	DO l=1,llm
	1034	!CDIR ON_ADB(index_i)
	1035	!CDIR ON_ADB(index_j)
	1036	!cdir NODEP
	1037	DO ig0=kstart,kend
	1038	i=index_i(ig0)
	1039	j=index_j(ig0)
	1040	pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
	1041	if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
	1042	ENDDO
	1043
	1044	IF (is_north_pole) then
	1045	DO i=1,iip1
	1046	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
	1047	ENDDO
	1048	ENDIF
	1049
	1050	IF (is_south_pole) then
	1051	DO i=1,iip1
	1052	pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
	1053	ENDDO
	1054	ENDIF
	1055
	1056	ENDDO
	1057	c$OMP END DO NOWAIT
	1058	ENDDO
	1059
	1060	c 65. champ u:
	1061	c ------------
	1062	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1063	DO l=1,llm
	1064	!CDIR ON_ADB(index_i)
	1065	!CDIR ON_ADB(index_j)
	1066	!cdir NODEP
	1067	do ig0=kstart,kend
	1068	i=index_i(ig0)
	1069	j=index_j(ig0)
	1070
	1071	if (i/=iim) then
	1072	pdufi(i,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
	1073	endif
	1074
	1075	if (i==1) then
	1076	pdufi(iim,j,l)=0.5*( zdufi2(ig0,l)
	1077	$ + zdufi2(ig0+iim-1,l))*cu(iim,j)
	1078	pdufi(iip1,j,l)=0.5(zdufi2(ig0,l)+zdufi2(ig0+1,l))cu(i,j)
	1079	endif
	1080
	1081	enddo
	1082
	1083	if (is_north_pole) then
	1084	DO i=1,iip1
	1085	pdufi(i,1,l) = 0.
	1086	ENDDO
	1087	endif
	1088
	1089	if (is_south_pole) then
	1090	DO i=1,iip1
	1091	pdufi(i,jjp1,l) = 0.
	1092	ENDDO
	1093	endif
	1094
	1095	ENDDO
	1096	c$OMP END DO NOWAIT
	1097
	1098	c 67. champ v:
	1099	c ------------
	1100
	1101	kstart=1
	1102	kend=klon
	1103
	1104	if (is_north_pole) kstart=2
	1105	if (is_south_pole) kend=klon-1-iim
	1106
	1107	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1108	DO l=1,llm
	1109	!CDIR ON_ADB(index_i)
	1110	!CDIR ON_ADB(index_j)
	1111	!cdir NODEP
	1112	do ig0=kstart,kend
	1113	i=index_i(ig0)
	1114	j=index_j(ig0)
	1115	pdvfi(i,j,l)=0.5(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))cv(i,j)
	1116	if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
	1117	$ zdvfi2(ig0+iim,l))
	1118	$ *cv(i,j)
	1119	enddo
	1120
	1121	ENDDO
	1122	c$OMP END DO NOWAIT
	1123
	1124
	1125	c 68. champ v pres des poles:
	1126	c ---------------------------
	1127	c v = U * cos(long) + V * SIN(long)
	1128
	1129	if (is_north_pole) then
	1130
	1131	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1132	DO l=1,llm
	1133
	1134	DO i=1,iim
	1135	pdvfi(i,1,l)=
	1136	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
	1137
	1138	pdvfi(i,1,l)=
	1139	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
	1140	ENDDO
	1141
	1142	pdvfi(iip1,1,l) = pdvfi(1,1,l)
	1143
	1144	ENDDO
	1145	c$OMP END DO NOWAIT
	1146
	1147	endif
	1148
	1149	if (is_south_pole) then
	1150
	1151	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1152	DO l=1,llm
	1153
	1154	DO i=1,iim
	1155	pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
	1156	$ +zdvfi(klon,l)*SIN(rlonv(i))
	1157
	1158	pdvfi(i,jjm,l)=
	1159	$ 0.5(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))cv(i,jjm)
	1160	ENDDO
	1161
	1162	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
	1163
	1164	ENDDO
	1165	c$OMP END DO NOWAIT
	1166
	1167	endif
	1168	c-----------------------------------------------------------------------
	1169
	1170	700 CONTINUE
	1171
	1172	firstcal = .FALSE.
	1173
	1174	#else
	1175	write(lunout,*)
	1176	& "calfis_p: for now can only work with parallel physics"
	1177	stop
	1178	#endif
[8]	1179	! of #ifdef CPP_PHYS
[1]	1180	RETURN
	1181	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format