Context Navigation

surf_heat_transp_mod.F90 @ 2176

Last change on this file since 2176 was 2111, checked in by emillour, 6 years ago

Generic model:

some fixes for the slab ocean. Still need to make it work in parallel.

Property svn:executable set to *

File size: 15.7 KB

Rev	Line
[1298]	1	!
	2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ocean_slab_mod.F90,v 1.3 2008-02-04 16:24:28 fairhead Exp $
	3	!
	4	MODULE surf_heat_transp_mod
	5
[1529]	6	IMPLICIT NONE
[1298]	7
[1529]	8	! Variables copied over from dyn3d dynamics:
	9	REAL,SAVE,ALLOCATABLE :: fext(:)
	10	REAL,SAVE,ALLOCATABLE :: unsairez(:)
	11	REAL,SAVE,ALLOCATABLE :: unsaire(:)
	12	REAL,SAVE,ALLOCATABLE :: cu(:)
	13	REAL,SAVE,ALLOCATABLE :: cv(:)
	14	REAL,SAVE,ALLOCATABLE :: cuvsurcv(:)
	15	REAL,SAVE,ALLOCATABLE :: cvusurcu(:)
	16	REAL,SAVE,ALLOCATABLE :: aire(:)
	17	REAL,SAVE :: apoln
	18	REAL,SAVE :: apols
	19	REAL,SAVE,ALLOCATABLE :: aireu(:)
	20	REAL,SAVE,ALLOCATABLE :: airev(:)
	21
	22	LOGICAL,SAVE :: alpha_var
	23	LOGICAL,SAVE :: slab_upstream
	24	REAL,SAVE,ALLOCATABLE :: zmasqu(:)
	25	REAL,SAVE,ALLOCATABLE :: zmasqv(:)
	26	REAL,SAVE,ALLOCATABLE :: unsfv(:)
	27	REAL,SAVE,ALLOCATABLE :: unsev(:)
	28	REAL,SAVE,ALLOCATABLE :: unsfu(:)
	29	REAL,SAVE,ALLOCATABLE :: unseu(:)
	30
	31	! Routines usable only by routines within this module:
	32	PRIVATE :: gr_fi_dyn, gr_dyn_fi
	33	! Routines from dyn3d, valid on global dynamics grid only:
[2111]	34	PRIVATE :: grad,diverg,gr_v_scal,gr_scal_v,gr_scal_u
[1529]	35
[1298]	36	CONTAINS
[1529]	37
	38	SUBROUTINE ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,&
	39	cu_,cuvsurcv_,cv_,cvusurcu_, &
	40	aire_,apoln_,apols_, &
	41	aireu_,airev_)
	42	USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
	43	IMPLICIT NONE
	44	! Transfer some variables from dyn3d dynamics
	45	INTEGER,INTENT(IN) :: ip1jm
	46	INTEGER,INTENT(IN) :: ip1jmp1
	47	REAL,INTENT(IN) :: unsairez_(ip1jm)
	48	REAL,INTENT(IN) :: fext_(ip1jm)
	49	REAL,INTENT(IN) :: unsaire_(ip1jmp1)
	50	REAL,INTENT(IN) :: cu_(ip1jmp1)
	51	REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
	52	REAL,INTENT(IN) :: cv_(ip1jm)
	53	REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
	54	REAL,INTENT(IN) :: aire_(ip1jmp1)
	55	REAL,INTENT(IN) :: apoln_
	56	REAL,INTENT(IN) :: apols_
	57	REAL,INTENT(IN) :: aireu_(ip1jmp1)
	58	REAL,INTENT(IN) :: airev_(ip1jm)
	59
	60	! Sanity check
	61	if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
	62	(ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
	63	write(,) "ini_surf_heat_transp Error: wrong array sizes"
	64	stop
	65	endif
	66
	67	allocate(unsairez(ip1jm))
	68	unsairez(:)=unsairez_(:)
	69	allocate(fext(ip1jm))
	70	fext(:)=fext_(:)
	71	allocate(unsaire(ip1jmp1))
	72	unsaire(:)=unsaire_(:)
	73	allocate(cu(ip1jmp1))
	74	cu(:)=cu_(:)
	75	allocate(cuvsurcv(ip1jm))
	76	cuvsurcv(:)=cuvsurcv_(:)
	77	allocate(cv(ip1jm))
	78	cv(:)=cv_(:)
	79	allocate(cvusurcu(ip1jmp1))
	80	cvusurcu(:)=cvusurcu_(:)
	81	allocate(aire(ip1jmp1))
	82	aire(:)=aire_(:)
	83	apoln=apoln_
	84	apols=apols_
	85	allocate(aireu(ip1jmp1))
	86	aireu(:)=aireu_(:)
	87	allocate(airev(ip1jm))
	88	airev(:)=airev_(:)
	89
	90	END SUBROUTINE ini_surf_heat_transp
	91
	92	SUBROUTINE ini_surf_heat_transp_mod
	93	USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
	94	IMPLICIT NONE
	95	INTEGER :: ip1jm, ip1jmp1
	96
	97	ip1jm=(nbp_lon+1)*(nbp_lat-1)
	98	ip1jmp1=(nbp_lon+1)*nbp_lat
	99
	100	allocate(zmasqu(ip1jmp1))
	101	allocate(zmasqv(ip1jm))
	102	allocate(unsfv(ip1jm))
	103	allocate(unsev(ip1jm))
	104	allocate(unsfu(ip1jmp1))
	105	allocate(unseu(ip1jmp1))
[1298]	106
[1529]	107	END SUBROUTINE ini_surf_heat_transp_mod
[1298]	108
[1529]	109	SUBROUTINE divgrad_phy(ngrid,nlevs,temp,delta)
[1298]	110
[1529]	111	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
[1298]	112	IMPLICIT NONE
	113
[1308]	114	INTEGER,INTENT(IN) :: ngrid, nlevs
	115	REAL,INTENT(IN) :: temp(ngrid,nlevs)
	116	REAL,INTENT(OUT) :: delta(ngrid,nlevs)
[1529]	117	REAL delta_2d((nbp_lon+1)*nbp_lat,nlevs)
[1308]	118	INTEGER :: ll
[1529]	119	REAL ghx((nbp_lon+1)*nbp_lat,nlevs)
	120	REAL ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
	121	INTEGER :: iip1,jjp1
	122
	123	iip1=nbp_lon+1
	124	jjp1=nbp_lat
[1298]	125
[1308]	126	CALL gr_fi_dyn(nlevs,ngrid,iip1,jjp1,temp,delta_2d)
[1298]	127	CALL grad(nlevs,delta_2d,ghx,ghy)
	128	DO ll=1,nlevs
	129	ghx(:,ll)=ghx(:,ll)*zmasqu
	130	! pas de diffusion zonale
	131	! ghx(:,ll)=0.
	132	ghy(:,ll)=ghy(:,ll)*zmasqv
	133	END DO
	134	CALL diverg(nlevs,ghx,ghy,delta_2d)
[1308]	135	CALL gr_dyn_fi(nlevs,iip1,jjp1,ngrid,delta_2d,delta)
[1298]	136
	137
	138	END SUBROUTINE divgrad_phy
	139
	140
	141
[1529]	142	SUBROUTINE init_masquv(ngrid,zmasq)
[1397]	143
[1529]	144	USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
[1298]	145	IMPLICIT NONE
	146
	147
[1308]	148	INTEGER,INTENT(IN) :: ngrid
	149	REAL zmasq(ngrid)
[1529]	150	REAL zmasq_2d((nbp_lon+1)*nbp_lat)
	151	REAL ff((nbp_lon+1)*(nbp_lat-1))
[1298]	152	REAL eps
	153	INTEGER i
[1529]	154	INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
	155
	156	iim=nbp_lon
	157	iip1=nbp_lon+1
	158	jjp1=nbp_lat
	159	ip1jm=(nbp_lon+1)*(nbp_lat-1)
	160	ip1jmp1=(nbp_lon+1)*nbp_lat
[1298]	161
	162	! Masques u,v
	163	zmasqu=1.
	164	zmasqv=1.
	165
[1308]	166	CALL gr_fi_dyn(1,ngrid,iip1,jjp1,zmasq,zmasq_2d)
[1298]	167
	168	DO i=1,ip1jmp1-1
	169	IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN
	170	zmasqu(i)=0.
	171	ENDIF
	172	END DO
	173	DO i=iip1,ip1jmp1,iip1
	174	zmasqu(i)=zmasqu(i-iim)
	175	END DO
	176	DO i=1,ip1jm
	177	IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN
	178	zmasqv(i)=0.
	179	END IF
	180	END DO
	181
	182
	183	! Coriolis (pour Ekman transp.)
	184	eps=1e-5
	185	! CALL getin('slab_eps',eps)
	186	! print *,'epsilon=',eps
	187	ff=fext*unsairez
	188	DO i=1,ip1jm
	189	unsev(i)=eps/(ff(i)ff(i)+eps*2)
	190	unsfv(i)=ff(i)/(ff(i)ff(i)+eps*2)
	191	ENDDO
	192	CALL gr_v_scal(1,unsfv,unsfu)
	193	CALL gr_v_scal(1,unsev,unseu)
	194	! Alpha variable?
	195	! alpha_var=.FALSE.
	196	! CALL getin('slab_alphav',alpha_var)
	197
	198
	199
	200	END SUBROUTINE init_masquv
	201
	202
	203
[1529]	204	SUBROUTINE slab_ekman2(ngrid,tx_phy,ty_phy,ts_phy,dt_phy)
[1298]	205
[1529]	206	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
	207	USE slab_ice_h, ONLY: noceanmx
[1298]	208
	209	IMPLICIT NONE
	210
[1308]	211	INTEGER,INTENT(IN) :: ngrid
[1298]	212	INTEGER ij
[1529]	213	REAL txv((nbp_lon+1)(nbp_lat-1)),fluxm((nbp_lon+1)(nbp_lat-1))
	214	REAL tyv((nbp_lon+1)*(nbp_lat-1))
	215	REAL fluxtm((nbp_lon+1)*(nbp_lat-1),noceanmx)
	216	REAL fluxtz((nbp_lon+1)*nbp_lat,noceanmx)
	217	REAL tyu((nbp_lon+1)nbp_lat),txu((nbp_lon+1)nbp_lat)
	218	REAL fluxz((nbp_lon+1)nbp_lat),fluxv((nbp_lon+1)nbp_lat)
	219	REAL dt((nbp_lon+1)nbp_lat,noceanmx),ts((nbp_lon+1)nbp_lat,noceanmx)
[1308]	220	REAL tx_phy(ngrid),ty_phy(ngrid)
	221	REAL dt_phy(ngrid,noceanmx),ts_phy(ngrid,noceanmx)
[1298]	222
[1529]	223	INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
[1298]	224
[1529]	225	iim=nbp_lon
	226	iip1=nbp_lon+1
	227	iip2=nbp_lon+2
[1537]	228	jjp1=nbp_lat
[1529]	229	ip1jm=(nbp_lon+1)(nbp_lat-1) ! = iip1jjm
	230	ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
	231	ip1jmp1=(nbp_lon+1)nbp_lat ! = iip1jjp1
[1298]	232
[2111]	233	! Convert taux,y to 2D scalar grid
	234	! north and south poles tx,ty no meaning
	235	tx_phy(1)=0.
	236	tx_phy(ngrid)=0.
	237	ty_phy(1)=0.
	238	ty_phy(ngrid)=0.
[1308]	239	CALL gr_fi_dyn(1,ngrid,iip1,jjp1,tx_phy,txu)
	240	CALL gr_fi_dyn(1,ngrid,iip1,jjp1,ty_phy,tyu)
[2111]	241
	242	! Divide taux,y by f or eps, and convert to 2D u,v grids
	243	! (Arakawa C grid)
	244	CALL gr_scal_v(1,txu,txv) ! wind stress at v points
	245	CALL gr_scal_v(1,tyu,tyv)
[1298]	246	fluxm=tyvunsev-txvunsfv
	247	! fluxm=-txv*unsfv
[2111]	248	! Zonal flux
	249	CALL gr_scal_u(1,txu,txu) ! wind stress at u points
	250	CALL gr_scal_u(1,tyu,tyu)
[1298]	251	fluxz=tyuunsfu+txuunseu
	252	! fluxz=tyu*unsfu
	253
[2111]	254	! Convert temperature to 2D grid
[1308]	255	CALL gr_fi_dyn(2,ngrid,iip1,jjp1,ts_phy,ts)
[1298]	256
	257	! Flux de masse
	258	fluxm=fluxmcvcuvsurcv*zmasqv
	259	fluxz=fluxzcucvusurcu*zmasqu
	260	! Flux de masse vertical
	261	DO ij=iip2,ip1jm
	262	fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
	263	ENDDO
	264	DO ij=iip1,ip1jmi1,iip1
	265	fluxv(ij+1)=fluxv(ij+iip1)
	266	END DO
	267	! Poles
	268	fluxv(1)=-SUM(fluxm(1:iim))
	269	fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
	270	fluxv=fluxv
	271
[2111]	272	! Meridional heat fluxes
[1298]	273	DO ij=1,ip1jm
[2111]	274	! centered scheme
[1298]	275	fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
	276	fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
	277	END DO
	278
[2111]	279	! Zonal heat fluxes
	280	! Schema upstream
	281	fluxtz(1:iip1,:)=0 ! no zonal heat flux at north pole
[1298]	282	DO ij=iip2,ip1jm
	283	IF (fluxz(ij).GE.0.) THEN
	284	fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
	285	fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
	286	ELSE
	287	fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
	288	fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
	289	ENDIF
	290	END DO
	291	DO ij=iip1*2,ip1jmp1,iip1
	292	fluxtz(ij,:)=fluxtz(ij-iim,:)
	293	END DO
	294
	295	! Calcul de dT
	296	DO ij=iip2,ip1jm
	297	dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
	298	+fluxtm(ij,:)-fluxtm(ij-iip1,:)
	299	IF (fluxv(ij).GT.0.) THEN
	300	dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2)
	301	dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2)
	302	ELSE
	303	dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1)
	304	dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1)
	305	ENDIF
	306	dt(ij,:)=dt(ij,:)/aire(ij)
	307	END DO
	308	DO ij=iip1,ip1jmi1,iip1
	309	dt(ij+1,:)=dt(ij+iip1,:)
	310	END DO
	311	! Pôles
	312	dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
	313	IF (fluxv(1).GT.0.) THEN
	314	dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2)
	315	dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2)
	316	ELSE
	317	dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1)
	318	dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1)
	319	ENDIF
	320	dt(1,:)=dt(1,:)/apoln
	321	dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
	322	IF (fluxv(ip1jmp1).GT.0.) THEN
	323	dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2)
	324	dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2)
	325	ELSE
	326	dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1)
	327	dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1)
	328	ENDIF
	329	dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols
	330	dt(2:iip1,1)=dt(1,1)
	331	dt(2:iip1,2)=dt(1,2)
	332	dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1)
	333	dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2)
	334
	335	! Retour grille physique
[1308]	336	CALL gr_dyn_fi(2,iip1,jjp1,ngrid,dt,dt_phy)
[1298]	337
	338
	339	END SUBROUTINE slab_ekman2
	340
[1529]	341	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	342
	343	SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
	344	! Transfer a variable on global "physics" grid to global "dynamics" grid
	345	IMPLICIT NONE
	346
	347	INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
	348	REAL,INTENT(IN) :: pfi(ngrid,nfield)
	349	REAL,INTENT(OUT) :: pdyn(im,jm,nfield)
	350
	351	INTEGER :: i,j,ifield,ig
	352
	353	DO ifield=1,nfield
	354	! Handle poles
	355	DO i=1,im
	356	pdyn(i,1,ifield)=pfi(1,ifield)
	357	pdyn(i,jm,ifield)=pfi(ngrid,ifield)
	358	ENDDO
	359	! Other points
	360	DO j=2,jm-1
	361	ig=2+(j-2)*(im-1)
	362	CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
	363	pdyn(im,j,ifield)=pdyn(1,j,ifield)
	364	ENDDO
	365	ENDDO ! of DO ifield=1,nfield
	366
	367	END SUBROUTINE gr_fi_dyn
	368
	369
	370
	371	SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
	372	! Transfer a variable on global "dynamics" grid to global "physics" grid
	373	IMPLICIT NONE
	374
	375	INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
	376	REAL,INTENT(IN) :: pdyn(im,jm,nfield)
	377	REAL,INTENT(OUT) :: pfi(ngrid,nfield)
	378
	379	INTEGER j,ifield,ig
	380
	381	! Sanity check:
	382	IF(ngrid.NE.2+(jm-2)*(im-1)) THEN
	383	WRITE(,) "gr_dyn_fi error, wrong sizes"
	384	STOP
	385	ENDIF
	386
	387	! Handle poles
	388	CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
	389	CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
	390	! Other points
	391	DO ifield=1,nfield
	392	DO j=2,jm-1
	393	ig=2+(j-2)*(im-1)
	394	CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
	395	ENDDO
	396	ENDDO
	397
	398	END SUBROUTINE gr_dyn_fi
	399
	400
	401
	402	SUBROUTINE grad(klevel,pg,pgx,pgy)
	403	! compute the covariant components x,y of the gradient of pg
	404	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
	405	IMPLICIT NONE
	406
	407	INTEGER,INTENT(IN) :: klevel
	408	REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
	409	REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
	410	REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
	411
	412	INTEGER :: l,ij
	413	INTEGER :: iim,iip1,ip1jm,ip1jmp1
	414
	415	iim=nbp_lon
	416	iip1=nbp_lon+1
	417	ip1jm=(nbp_lon+1)(nbp_lat-1) ! = iip1jjm
	418	ip1jmp1=(nbp_lon+1)nbp_lat ! = iip1jjp1
	419
	420	DO l=1,klevel
	421	DO ij=1,ip1jmp1-1
	422	pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
	423	ENDDO
	424	! correction for pgx(ip1,j,l) ...
	425	! ... pgx(iip1,j,l)=pgx(1,j,l) ...
	426	DO ij=iip1,ip1jmp1,iip1
	427	pgx(ij,l)=pgx(ij-iim,l)
	428	ENDDO
	429	DO ij=1,ip1jm
	430	pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
	431	ENDDO
	432	ENDDO
	433
	434	END SUBROUTINE grad
	435
	436
	437
	438	SUBROUTINE diverg(klevel,x,y,div)
	439	! compute the divergence of a vector field of components
	440	! x,y. y and y being covriant components
	441	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
	442	IMPLICIT NONE
	443
	444	INTEGER,INTENT(IN) :: klevel
	445	REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
	446	REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
	447	REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
	448
	449	INTEGER :: l,ij
	450	INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
	451
	452	REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
	453	REAL :: sumypn,sumyps
	454	REAL,EXTERNAL :: SSUM
	455
	456	iim=nbp_lon
	457	iip1=nbp_lon+1
	458	iip2=nbp_lon+2
	459	ip1jm=(nbp_lon+1)(nbp_lat-1) ! = iip1jjm
	460	ip1jmp1=(nbp_lon+1)nbp_lat ! = iip1jjp1
	461	ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
	462
	463	DO l=1,klevel
	464	DO ij=iip2,ip1jm-1
	465	div(ij+1,l)= &
	466	cvusurcu(ij+1)x(ij+1,l)-cvusurcu(ij)x(ij,l)+ &
	467	cuvsurcv(ij-iim)y(ij-iim,l)-cuvsurcv(ij+1)y(ij+1,l)
	468	ENDDO
	469	! correction for div(1,j,l) ...
	470	! ... div(1,j,l)= div(iip1,j,l) ...
	471	DO ij=iip2,ip1jm,iip1
	472	div(ij,l)=div(ij+iim,l)
	473	ENDDO
	474	! at the poles
	475	DO ij=1,iim
	476	aiy1(ij)=cuvsurcv(ij)*y(ij,l)
	477	aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
	478	ENDDO
	479	sumypn=SSUM(iim,aiy1,1)/apoln
	480	sumyps=SSUM(iim,aiy2,1)/apols
	481	DO ij=1,iip1
	482	div(ij,l)=-sumypn
	483	div(ij+ip1jm,l)=sumyps
	484	ENDDO
	485	ENDDO ! of DO l=1,klevel
	486
	487	!!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
	488
	489	DO l=1,klevel
	490	DO ij=iip2,ip1jm
	491	div(ij,l)=div(ij,l)*unsaire(ij)
	492	ENDDO
	493	ENDDO
	494
	495	END SUBROUTINE diverg
	496
	497
	498
	499	SUBROUTINE gr_v_scal(nx,x_v,x_scal)
	500	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
	501	IMPLICIT NONE
	502
	503	INTEGER,INTENT(IN) :: nx
	504	REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
	505	REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
	506
	507	INTEGER :: l,ij
	508	INTEGER :: iip1,iip2,ip1jm,ip1jmp1
	509
	510	iip1=nbp_lon+1
	511	iip2=nbp_lon+2
	512	ip1jm=(nbp_lon+1)(nbp_lat-1) ! = iip1jjm
	513	ip1jmp1=(nbp_lon+1)nbp_lat ! = iip1jjp1
	514
	515	DO l=1,nx
	516	DO ij=iip2,ip1jm
	517	x_scal(ij,l)= &
	518	(airev(ij-iip1)x_v(ij-iip1,l)+airev(ij)x_v(ij,l)) &
	519	/(airev(ij-iip1)+airev(ij))
	520	ENDDO
	521	DO ij=1,iip1
	522	x_scal(ij,l)=0.
	523	ENDDO
	524	DO ij=ip1jm+1,ip1jmp1
	525	x_scal(ij,l)=0.
	526	ENDDO
	527	ENDDO
	528
	529	END SUBROUTINE gr_v_scal
[2111]	530
	531	SUBROUTINE gr_scal_v(nx,x_scal,x_v)
	532	! convert values from scalar points to v points on C-grid
	533	! used to compute wind stress at V points
	534	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
	535	IMPLICIT NONE
	536
	537	INTEGER,INTENT(IN) :: nx ! number of levels or fields
	538	REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
	539	REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
	540
	541	INTEGER :: l,ij
	542	INTEGER :: iip1,ip1jm
	543
	544	iip1=nbp_lon+1
	545	ip1jm=(nbp_lon+1)(nbp_lat-1) ! = iip1jjm
	546
	547	DO l=1,nx
	548	DO ij=1,ip1jm
	549	x_v(ij,l)= &
	550	(cu(ij)cvusurcu(ij)x_scal(ij,l)+ &
	551	cu(ij+iip1)cvusurcu(ij+iip1)x_scal(ij+iip1,l)) &
	552	/(cu(ij)cvusurcu(ij)+cu(ij+iip1)cvusurcu(ij+iip1))
	553	ENDDO
	554	ENDDO
	555
	556	END SUBROUTINE gr_scal_v
	557
	558	SUBROUTINE gr_scal_u(nx,x_scal,x_u)
	559	! convert values from scalar points to U points on C-grid
	560	! used to compute wind stress at U points
	561	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
	562	IMPLICIT NONE
	563
	564	INTEGER,INTENT(IN) :: nx
	565	REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx)
	566	REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
	567
	568	INTEGER :: l,ij
	569	INTEGER :: iip1,jjp1,ip1jmp1
	570
	571	iip1=nbp_lon+1
	572	jjp1=nbp_lat
	573	ip1jmp1=(nbp_lon+1)nbp_lat ! = iip1jjp1
	574
	575	DO l=1,nx
	576	DO ij=1,ip1jmp1-1
	577	x_u(ij,l)= &
	578	(aire(ij)x_scal(ij,l)+aire(ij+1)x_scal(ij+1,l)) &
	579	/(aire(ij)+aire(ij+1))
	580	ENDDO
	581	ENDDO
	582
	583	CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1)
	584
	585	END SUBROUTINE gr_scal_u
[1529]	586
[1298]	587	END MODULE surf_heat_transp_mod
	588
	589
	590

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

Download in other formats:

Original Format