Index: DZ4/branches/LMDZ4-dev/libf/dyn3dpar/conf_guide.F =================================================================== --- /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/conf_guide.F (revision 1172) +++ (revision ) @@ -1,57 +1,0 @@ -! -! $Header$ -! -c -c - SUBROUTINE conf_guide -c - use IOIPSL - use getparam - IMPLICIT NONE - -c----------------------------------------------------------------------- -c Parametres de controle du run: -c----------------------------------------------------------------------- -#include "guide.h" - - call getpar('guide.eff') - - call getpar('online',1,online,'Index de controle du guide') - CALL getpar('ncep',.false.,ncep,'Coordonnee vert NCEP ou ECMWF') - CALL getpar('guide_modele',.false.,guide_modele, - $ 'guidage sur niveaux modele') - CALL getpar('ini_anal',.false.,ini_anal,'Initial = analyse') - CALL getpar('ok_invertp',.true.,ok_invertp,'niveaux p inverses') - - CALL getpar('guide_u',.true.,guide_u,'guidage de u') - CALL getpar('guide_v',.true.,guide_v,'guidage de v') - CALL getpar('guide_T',.true.,guide_T,'guidage de T') - CALL getpar('guide_P',.true.,guide_P,'guidage de P') - CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q') - CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R') - CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim') - -c Constantes de rappel. Unite : fraction de jour - CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u') - CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u') - CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v') - CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v') - CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T') - CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T') - CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q') - CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q') - CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P') - CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P') - CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie') - -c Latitude min et max pour le rappel. -c dans le cas ou on 'a les analyses que sur une bande de latitudes. - CALL getpar('lat_min_guide',-90.,lat_min_guide - s ,'Latitude minimum pour le guidage ') - CALL getpar('lat_max_guide', 90.,lat_max_guide - s ,'Latitude maximum pour le guidage ') - - - CALL getpar - - end Index: DZ4/branches/LMDZ4-dev/libf/dyn3dpar/guide.h =================================================================== --- /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/guide.h (revision 1172) +++ (revision ) @@ -1,42 +1,0 @@ -! -! $Header$ -! - real tau_min_u,tau_max_u - real tau_min_v,tau_max_v - real tau_min_T,tau_max_T - real tau_min_q,tau_max_q - real tau_min_P,tau_max_P - real aire_min,aire_max - - - logical guide_u,guide_v,guide_T,guide_Q,guide_P,guide_modele - logical guide_BL,guide_hr,gamma4 - - real lat_min_guide,lat_max_guide - - -c data tau_min_u,tau_max_u/0.02,10./ -c data tau_min_v,tau_max_v/0.02,10./ -c data tau_min_T,tau_max_T/0.02,10./ -c data tau_min_q,tau_max_q/0.02,10./ -c data tau_min_P,tau_max_P/0.02,10./ -c - LOGICAL ncep,ini_anal,ok_invertp - integer online - -c data online/1/ -c data ncep,ini_anal/.false.,.true./ - - common/comguide/ - s tau_min_u,tau_max_u, - s tau_min_v,tau_max_v, - s tau_min_T,tau_max_T, - s tau_min_q,tau_max_q, - s tau_min_P,tau_max_P, - s aire_min,aire_max, - s lat_min_guide,lat_max_guide, - s ncep,ini_anal, - s online, - s guide_u,guide_v,guide_T,guide_Q,guide_P, - s guide_modele,ok_invertp, - s guide_hr,guide_BL,gamma4 Index: DZ4/branches/LMDZ4-dev/libf/dyn3dpar/guide_p.F =================================================================== --- /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/guide_p.F (revision 1172) +++ (revision ) @@ -1,590 +1,0 @@ -! -! $Id$ -! - subroutine guide_pp(itau,ucov,vcov,teta,q,masse,ps) - USE parallel - use netcdf - - IMPLICIT NONE - -c ...... Version du 10/01/98 .......... - -c avec coordonnees verticales hybrides -c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 ) - -c======================================================================= -c -c Auteur: F.Hourdin -c ------- -c -c Objet: -c ------ -c -c GCM LMD nouvelle grille -c -c======================================================================= - -c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv -c et possibilite d'appeler une fonction f(y) a derivee tangente -c hyperbolique a la place de la fonction a derivee sinusoidale. - -c ... Possibilite de choisir le shema de Van-leer pour l'advection de -c q , en faisant iadv = 10 dans traceur (29/04/97) . -c -c----------------------------------------------------------------------- -c Declarations: -c ------------- - -#include "dimensions.h" -#include "paramet.h" -#include "comconst.h" -#include "comdissnew.h" -#include "comvert.h" -#include "comgeom.h" -#include "logic.h" -#include "temps.h" -#include "control.h" -#include "ener.h" -#include "netcdf.inc" -#include "description.h" -#include "serre.h" -#include "tracstoke.h" -#include "guide.h" - - -c variables dynamiques - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants - REAL teta(ip1jmp1,llm) ! temperature potentielle - REAL q(ip1jmp1,llm) ! temperature potentielle - REAL ps(ip1jmp1) ! pression au sol - REAL masse(ip1jmp1,llm) ! masse d'air - -c common passe pour des sorties - real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm) - common/comdxdy/dxdys,dxdyu,dxdyv - -c variables dynamiques pour les reanalyses. - REAL ucovrea1(ip1jmp1,llm),vcovrea1(ip1jm,llm) !vts cov reas - REAL tetarea1(ip1jmp1,llm) ! temp pot reales - REAL qrea1(ip1jmp1,llm) ! temp pot reales - REAL masserea1(ip1jmp1,llm) ! masse - REAL psrea1(ip1jmp1) ! ps - REAL ucovrea2(ip1jmp1,llm),vcovrea2(ip1jm,llm) !vts cov reas - REAL tetarea2(ip1jmp1,llm) ! temp pot reales - REAL qrea2(ip1jmp1,llm) ! temp pot reales - REAL masserea2(ip1jmp1,llm) ! masse - REAL psrea2(ip1jmp1) ! ps - real latmin - - real alpha_q(ip1jmp1) - real alpha_T(ip1jmp1),alpha_P(ip1jmp1) - real alpha_u(ip1jmp1),alpha_v(ip1jm) - real alpha_pcor(llm) - real dday_step,toto,reste,itau_test - INTEGER step_rea,count_no_rea - -c real aire_min,aire_max - integer ilon,ilat - real factt,ztau(ip1jmp1) - - INTEGER itau,ij,l,i,j - integer ncidpl,varidpl,nlev,status - integer rcod,rid - real ditau,tau,a - save nlev - -c TEST SUR QSAT - real p(ip1jmp1,llmp1),pk(ip1jmp1,llm),pks(ip1jmp1) - real pkf(ip1jmp1,llm) - real pres(ip1jmp1,llm) - REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) - - real qsat(ip1jmp1,llm) - real unskap - real tnat(ip1jmp1,llm) -ccccccccccccccccc - - - LOGICAL first - save first - data first/.true./ - - save ucovrea1,vcovrea1,tetarea1,masserea1,psrea1,qrea1 - save ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,qrea2 - - save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test - save alpha_pcor - save step_rea,count_no_rea - - character*10 file - integer igrads - real dtgrads - save igrads,dtgrads - data igrads,dtgrads/2,100./ - integer :: ijb,ije,jjn - -C----------------------------------------------------------------------- -c calcul de l'humidite saturante -C----------------------------------------------------------------------- - ijb=ij_begin - ije=ij_end - jjn=jj_nb - - CALL pression_p( ip1jmp1, ap, bp, ps, p ) - call massdair_p(p,masse) - CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf) - tnat(ijb:ije,:)=pk(ijb:ije,:)*teta(ijb:ije,:)/cpp - unskap = 1./ kappa - pres(ijb:ije,:)=preff*(pk(ijb:ije,:)/cpp)**unskap - call q_sat(iip1*jjn*llm,tnat(ijb:ije,:),pres(ijb:ije,:), - . qsat(ijb:ije,:)) - -C----------------------------------------------------------------------- - -c----------------------------------------------------------------------- -c initialisations pour la lecture des reanalyses. -c alpha determine la part des injections de donnees a chaque etape -c alpha=1 signifie pas d'injection -c alpha=0 signifie injection totale -c----------------------------------------------------------------------- - - if(online.eq.-1) then - return - endif - - if (first) then - - print*,'initialisation du guide ' - call conf_guide - print*,'apres conf_guide' - - file='guide' - - if (mpi_rank==0) then - if (ok_gradsfile) then - call inigrads(igrads,iip1 - s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi - s ,llm,presnivs,1. - s ,dtgrads,file,'dyn_zon ') - endif !ok_gradsfile - endif - - print* - s ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)' - - if(online.eq.-1) return - if (online.eq.1) then - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c Constantes de temps de rappel en jour -c 0.1 c'est en gros 2h30. -c 1e10 est une constante infinie donc en gros pas de guidage -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c coordonnees du centre du zoom - call coordij(clon,clat,ilon,ilat) -c aire de la maille au centre du zoom - aire_min=aire(ilon+(ilat-1)*iip1) -c aire maximale de la maille - aire_max=0. - do ij=1,ip1jmp1 - aire_max=max(aire_max,aire(ij)) - enddo -C factt = pas de temps en fraction de jour - factt=dtvr*iperiod/daysec - -c subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha) - call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v) - call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u) - call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T) - call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P) - call tau2alpha(1,iip1,jjp1,factt,tau_min_q,tau_max_q,alpha_q) - - call dump2d(iip1,jjp1,aire,'AIRE MAILLe ') - call dump2d(iip1,jjp1,alpha_u,'COEFF U ') - call dump2d(iip1,jjp1,alpha_T,'COEFF T ') - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c Cas ou on force exactement par les variables analysees - else - alpha_T=0. - alpha_u=0. - alpha_v=0. - alpha_P=0. -c physic=.false. - endif -c correction de rappel dans couche limite -c F.Codron - if (guide_BL) then - alpha_pcor(:)=1. - else - do l=1,llm - alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2. - enddo - endif - itau_test=1001 - step_rea=1 - count_no_rea=0 - ncidpl=-99 -c itau_test montre si l'importation a deja ete faite au rang itau -c lecture d'un fichier netcdf pour determiner le nombre de niveaux - IF (mpi_rank==0) THEN - if (guide_modele) then - if (ncidpl.eq.-99) rcod=nf90_open('apbp.nc',Nf90_NOWRITe, - $ ncidpl) - else - if (guide_u) then - if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe, - $ ncidpl) - endif -c - if (guide_v) then - if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite, - $ ncidpl) - endif -c - if (guide_T) then - if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite, - $ ncidpl) - endif -c - if (guide_Q) then - if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, - $ ncidpl) - endif -c - endif !guide_modele -cAJ-b endif !mpi_rank - if (ncep) then - status=NF_INQ_DIMID(ncidpl,'LEVEL',rid) - else - status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid) - endif - status=NF_INQ_DIMLEN(ncidpl,rid,nlev) - print *,'nlev guide', nlev - rcod = nf90_close(ncidpl) -cAJ-e - endif !mpi_rank -cAJ-e -c Lecture du premier etat des reanalyses. - call Gather_Field(ps,ip1jmp1,1,0) - - if (mpi_rank==0) then - - call read_reanalyse(1,ps - s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev) - qrea2(:,:)=max(qrea2(:,:),0.1) - - endif - - call Broadcast_Field(ucovrea2,ip1jmp1,llm,0) - call Broadcast_Field(vcovrea2,ip1jm,llm,0) - call Broadcast_Field(tetarea2,ip1jmp1,llm,0) - call Broadcast_Field(qrea2,ip1jmp1,llm,0) - call Broadcast_Field(masserea2,ip1jmp1,llm,0) - call Broadcast_Field(psrea2,ip1jmp1,1,0) - - - -c----------------------------------------------------------------------- -c Debut de l'integration temporelle: -c ---------------------------------- - - endif ! first -c -C----------------------------------------------------------------------- -C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS: -C----------------------------------------------------------------------- - - ditau=real(itau) - DDAY_step=real(day_step) - write(*,*)'ditau,dday_step' - write(*,*)ditau,dday_step - toto=4*ditau/dday_step - reste=toto-aint(toto) -c write(*,*)'toto,reste',toto,reste - - if (reste.eq.0.) then - if (itau_test.eq.itau) then - write(*,*)'deuxieme passage de advreel a itau=',itau - stop - else - vcovrea1(:,:)=vcovrea2(:,:) - ucovrea1(:,:)=ucovrea2(:,:) - tetarea1(:,:)=tetarea2(:,:) - qrea1(:,:)=qrea2(:,:) - - print*,'LECTURE REANALYSES, pas ',step_rea - s ,'apres ',count_no_rea,' non lectures' - step_rea=step_rea+1 - itau_test=itau - - call Gather_Field(ps,ip1jmp1,1,0) - - if (mpi_rank==0) then - call read_reanalyse(step_rea,ps - s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev) - qrea2(:,:)=max(qrea2(:,:),0.1) - endif - - call Broadcast_Field(ucovrea2,ip1jmp1,llm,0) - call Broadcast_Field(vcovrea2,ip1jm,llm,0) - call Broadcast_Field(tetarea2,ip1jmp1,llm,0) - call Broadcast_Field(qrea2,ip1jmp1,llm,0) - call Broadcast_Field(masserea2,ip1jmp1,llm,0) - call Broadcast_Field(psrea2,ip1jmp1,1,0) - - factt=dtvr*iperiod/daysec - ztau(:)=factt/max(alpha_T(:),1.e-10) - - if (mpi_rank==0) then - if (ok_gradsfile) then - call wrgrads(igrads,1,aire ,'aire ','aire ' ) - call wrgrads(igrads,1,dxdys ,'dxdy ','dxdy ' ) - call wrgrads(igrads,1,alpha_u,'au ','au ' ) - call wrgrads(igrads,1,alpha_T,'at ','at ' ) - call wrgrads(igrads,1,ztau,'taut ','taut ' ) - call wrgrads(igrads,llm,ucov,'u ','u ' ) - call wrgrads(igrads,llm,ucovrea2,'ua ','ua ' ) - call wrgrads(igrads,llm,teta,'T ','T ' ) - call wrgrads(igrads,llm,tetarea2,'Ta ','Ta ' ) - call wrgrads(igrads,llm,qrea2,'Qa ','Qa ' ) - call wrgrads(igrads,llm,q,'Q ','Q ' ) - call wrgrads(igrads,llm,qsat,'QSAT ','QSAT ' ) - endif !(ok_gradsfile) then - endif - - endif - else - count_no_rea=count_no_rea+1 - endif - -C----------------------------------------------------------------------- -c Guidage -c x_gcm = a * x_gcm + (1-a) * x_reanalyses -C----------------------------------------------------------------------- - - if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE' - - ditau=real(itau) - dday_step=real(day_step) - - - tau=4*ditau/dday_step - tau=tau-aint(tau) - -c ucov - ijb=ij_begin - ije=ij_end - - if (guide_u) then - do l=1,llm - do ij=ijb,ije - a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l) - ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij) - $ *alpha_pcor(l)*a - if (first.and.ini_anal) ucov(ij,l)=a - enddo - enddo - endif - -c teta - if (guide_T) then - do l=1,llm - do ij=ijb,ije - a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l) - teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij) - $ *alpha_pcor(l)*a - if (first.and.ini_anal) teta(ij,l)=a - enddo - enddo - endif - -c P - if (guide_P) then - do ij=ijb,ije - a=(1.-tau)*psrea1(ij)+tau*psrea2(ij) - ps(ij)=(1.-alpha_P(ij))*ps(ij)+alpha_P(ij)*a - if (first.and.ini_anal) ps(ij)=a - enddo - CALL pression_p(ip1jmp1,ap,bp,ps,p) - CALL massdair_p(p,masse) - endif - - -c q - if (guide_Q) then - do l=1,llm - do ij=ijb,ije - a=(1.-tau)*qrea1(ij,l)+tau*qrea2(ij,l) -c hum relative en % -> hum specif - a=qsat(ij,l)*a*0.01 - q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij) - $ *alpha_pcor(l)*a - if (first.and.ini_anal) q(ij,l)=a - enddo - enddo - endif - -c vcov - if (pole_sud) ije=ij_end-iip1 - - if (guide_v) then - do l=1,llm - do ij=ijb,ije - a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l) - vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij) - $ *alpha_pcor(l)*a - if (first.and.ini_anal) vcov(ij,l)=a - enddo - enddo - endif - -c call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ') -c call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ') -c call dump2d(iip1,jjp1,teta,'TETA ') - - first=.false. - - return - end - -c======================================================================= - subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha) -c======================================================================= - - implicit none - -#include "dimensions.h" -#include "paramet.h" -#include "comconst.h" -#include "comgeom2.h" -#include "guide.h" -#include "serre.h" - -c arguments : - integer type - integer pim,pjm - real factt,taumin,taumax,dxdymin,dxdymax - real dxdy_,alpha(pim,pjm) - real dxdy_min,dxdy_max - -c local : - real alphamin,alphamax,gamma,xi - save gamma - integer i,j,ilon,ilat - - logical first - save first - data first/.true./ - - real cus(iip1,jjp1),cvs(iip1,jjp1) - real cuv(iip1,jjm),cvu(iip1,jjp1) - real zdx(iip1,jjp1),zdy(iip1,jjp1) - - real zlat - real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm) - common/comdxdy/dxdys,dxdyu,dxdyv - - if (first) then - do j=2,jjm - do i=2,iip1 - zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j)) - enddo - zdx(1,j)=zdx(iip1,j) - enddo - do j=2,jjm - do i=1,iip1 - zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j)) - enddo - enddo - do i=1,iip1 - zdx(i,1)=zdx(i,2) - zdx(i,jjp1)=zdx(i,jjm) - zdy(i,1)=zdy(i,2) - zdy(i,jjp1)=zdy(i,jjm) - enddo - do j=1,jjp1 - do i=1,iip1 - dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j)) - enddo - enddo - do j=1,jjp1 - do i=1,iim - dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j)) - enddo - dxdyu(iip1,j)=dxdyu(1,j) - enddo - do j=1,jjm - do i=1,iip1 - dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j)) - enddo - enddo - -c call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ') -c call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ') - call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ') - -c coordonnees du centre du zoom - call coordij(clon,clat,ilon,ilat) -c aire de la maille au centre du zoom - dxdy_min=dxdys(ilon,ilat) -c dxdy maximale de la maille - dxdy_max=0. - do j=1,jjp1 - do i=1,iip1 - dxdy_max=max(dxdy_max,dxdys(i,j)) - enddo - enddo - - if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then - print*,'ATTENTION modele peu zoome' - print*,'ATTENTION on prend une constante de guidage cste' - gamma=0. - else - gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min) - print*,'gamma=',gamma - if (gamma.lt.1.e-5) then - print*,'gamma =',gamma,'<1e-5' - stop - endif - gamma=log(0.5)/log(gamma) - if (gamma4) then - gamma=min(gamma,4.) - endif - print*,'gamma=',gamma - endif - endif - - alphamin=factt/taumax - alphamax=factt/taumin - - do j=1,pjm - do i=1,pim - if (type.eq.1) then - dxdy_=dxdys(i,j) - zlat=rlatu(j)*180./pi - elseif (type.eq.2) then - dxdy_=dxdyu(i,j) - zlat=rlatu(j)*180./pi - elseif (type.eq.3) then - dxdy_=dxdyv(i,j) - zlat=rlatv(j)*180./pi - endif - if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then -c pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin - alpha(i,j)=alphamin - else - xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma - xi=min(xi,1.) - if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then - alpha(i,j)=xi*alphamin+(1.-xi)*alphamax - else - alpha(i,j)=0. - endif - endif - enddo - enddo - - - return - end Index: DZ4/branches/LMDZ4-dev/libf/dyn3dpar/read_reanalyse.F =================================================================== --- /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/read_reanalyse.F (revision 1172) +++ (revision ) @@ -1,580 +1,0 @@ -! -! $Header$ -! -c -c - subroutine read_reanalyse(timestep,psi - s ,u,v,t,q,masse,ps,mode,nlevnc) - -c mode=0 variables naturelles -c mode=1 variabels GCM - - USE parallel - use netcdf -c ----------------------------------------------------------------- -c Declarations -c ----------------------------------------------------------------- - IMPLICIT NONE - -c common -c ------ - -#include "netcdf.inc" -#include "dimensions.h" -#include "paramet.h" -#include "comgeom.h" -#include "comvert.h" -#include "guide.h" - -c arguments -c --------- - integer nlevnc - integer timestep,mode - - real psi(iip1,jjp1) - real u(iip1,jjp1,llm),v(iip1,jjm,llm) - real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm) - real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm) - -c local -c ----- - integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps - integer ncidpl - integer varidpl,ncidQ,varidQ,varidap,varidbp - save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps - save ncidpl,varidap,varidbp - save varidpl,ncidQ,varidQ - - real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc) - real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1) - real Qnc(iip1,jjp1,nlevnc) - real plnc(iip1,jjp1,nlevnc) - real apnc(nlevnc),bpnc(nlevnc) - - integer start(4),count(4),status - integer i,j,l - - real rcode - logical first - save first - INTEGER ierr - - data first/.true./ - -c ----------------------------------------------------------------- -c Initialisation de la lecture des fichiers -c ----------------------------------------------------------------- - if (first) then - ncidpl=-99 - print*,'Intitialisation de read reanalsye' - -c Niveaux de pression si non constants - if (guide_modele) then - print *,'Vous êtes entrain de lire des données sur - . niveaux modèle' - rcode=nf90_open('apbp.nc',nf90_nowrite,ncidpl) - rcode = nf90_inq_varid(ncidpl, 'AP', varidap) - rcode = nf90_inq_varid(ncidpl, 'BP', varidbp) - print*,'ncidpl,varidap',ncidpl,varidap - endif - -c Vent zonal - if (guide_u) then - rcode=nf90_open('u.nc',nf90_nowrite,ncidu) - rcode = nf90_inq_varid(ncidu, 'UWND', varidu) - print*,'ncidu,varidu',ncidu,varidu - if (ncidpl.eq.-99) ncidpl=ncidu - endif - -c Vent meridien - if (guide_v) then - rcode=nf90_open('v.nc',nf90_nowrite,ncidv) - rcode = nf90_inq_varid(ncidv, 'VWND', varidv) - print*,'ncidv,varidv',ncidv,varidv - if (ncidpl.eq.-99) ncidpl=ncidv - endif - -c Temperature - if (guide_T) then - rcode=nf90_open('T.nc',nf90_nowrite,ncidt) - rcode = nf90_inq_varid(ncidt, 'AIR', varidt) - print*,'ncidt,varidt',ncidt,varidt - if (ncidpl.eq.-99) ncidpl=ncidt - endif - -c Humidite - if (guide_Q) then - rcode=nf90_open('hur.nc',nf90_nowrite,ncidQ) - rcode = nf90_inq_varid(ncidQ, 'RH', varidQ) - print*,'ncidQ,varidQ',ncidQ,varidQ - if (ncidpl.eq.-99) ncidpl=ncidQ - endif - -c Pression de surface - if ((guide_P).OR.(guide_modele)) then - rcode=nf90_open('ps.nc',nf90_nowrite,ncidps) - rcode = nf90_inq_varid(ncidps, 'SP', varidps) - print*,'ncidps,varidps',ncidps,varidps - endif - -c Coordonnee verticale - if (.not.guide_modele) then - if (ncep) then - print*,'Vous etes entrain de lire des donnees NCEP' - rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl) - else - print*,'Vous etes entrain de lire des donnees ECMWF' - rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl) - endif - print*,'ncidpl,varidpl',ncidpl,varidpl - endif - print*,'ncidu,varidpl',ncidu,varidpl - - endif - print*,'ok1' - -c ----------------------------------------------------------------- -c lecture des champs u, v, T, ps -c ----------------------------------------------------------------- - -c dimensions pour les champs scalaires et le vent zonal -c ----------------------------------------------------- - - start(1)=1 - start(2)=1 - start(3)=1 - start(4)=timestep - - count(1)=iip1 - count(2)=jjp1 - count(3)=nlevnc - count(4)=1 - -c mise a zero des tableaux -c ------------------------ - unc(:,:,:)=0. - vnc(:,:,:)=0. - tnc(:,:,:)=0. - Qnc(:,:,:)=0. - plnc(:,:,:)=0. - -c Vent zonal -c ---------- - - if (guide_u) then - print*,'avant la lecture de UNCEP nd de niv:',nlevnc - -#ifdef NC_DOUBLE - status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc) -#else - status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc) -#endif - print*,'WARNING!!! Correction bidon pour palier a un ' - print*,'probleme dans la creation des fichiers nc' - call correctbid(iim,jjp1*nlevnc,unc) - call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ') - endif - -c Temperature -c ----------- - - print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start - print*,'count=',count - if (guide_T) then -#ifdef NC_DOUBLE - status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc) -#else - status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc) -#endif - call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 AAA ') - call correctbid(iim,jjp1*nlevnc,tnc) - call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ') - endif - -c Humidite -c -------- - - if (guide_Q) then -#ifdef NC_DOUBLE - status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,Qnc) -#else - status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc) -#endif - call correctbid(iim,jjp1*nlevnc,Qnc) - call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ') - endif - - count(2)=jjm -c Vent meridien -c ------------- - - if (guide_v) then -#ifdef NC_DOUBLE - status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc) -#else - status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc) -#endif - call correctbid(iim,jjm*nlevnc,vnc) - call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ') - endif - - start(3)=timestep - start(4)=0 - count(2)=jjp1 - count(3)=1 - count(4)=0 - -c Pression de surface -c ------------------- - psnc(:,:)=0. - if ((guide_P).OR.(guide_modele)) then -#ifdef NC_DOUBLE - status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc) -#else - status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc) -#endif - call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ') - call correctbid(iim,jjp1,psnc) - endif - -c Calcul de la pression aux differents niveaux -c -------------------------------------------- - if (guide_modele) then -#ifdef NC_DOUBLE - status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc) - status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc) -#else - status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc) - status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc) -#endif - else -#ifdef NC_DOUBLE - status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc) -#else - status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc) -#endif -c conversion en Pascals - apnc=apnc*100. - bpnc(:)=0. - endif - do i=1,iip1 - do j=1,jjp1 - do l=1,nlevnc - plnc(i,j,l)=apnc(l)+bpnc(l)*psnc(i,j) - enddo - enddo - enddo - if (first) then - print*,'Verification ordre niveaux verticaux' - print*,'LMDZ :' - do l=1,llm - print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. - $ +psi(1,jjp1)*(bp(l)+bp(l+1))/2. - enddo - print*,'Fichiers guidage' - do l=1,nlevnc - print*,'PL(',l,')=',plnc(1,1,l) - enddo - if (guide_u) then - do l=1,nlevnc - print*,'U(',l,')=',unc(1,1,l) - enddo - endif - if (guide_T) then - do l=1,nlevnc - print*,'T(',l,')=',tnc(1,1,l) - enddo - endif - print *,'inversion de l''ordre: ok_invertp=',ok_invertp - endif - -c ----------------------------------------------------------------- -c Interpollation verticale sur les niveaux modele -c ----------------------------------------------------------------- - call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,plnc,u,v,t,Q - s ,ps,masse,pk) - - call dump2d(iip1,jjm,v,'V COUCHE APRES ') - - -c ----------------------------------------------------------------- -c Passage aux variables du modele (vents covariants, temperature -c potentielle et humidite specifique) -c ----------------------------------------------------------------- - call nat2gcm(u,v,t,Q,pk,u,v,t,Q) - print*,'TIMESTEP ',timestep - if(mode.ne.1) stop'mode pas egal 0' -c call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ') - -c Lignes introduites a une epoque pour un probleme oublie... -c do l=1,llm -c do i=1,iip1 -c v(i,1,l)=0. -c v(i,jjm,l)=0. -c enddo -c enddo - first=.false. - - return - end - - -c=========================================================================== - subroutine reanalyse2nat(nlevnc,psi - s ,unc,vnc,tnc,qnc,psnc,plnc,u,v,t,q - s ,ps,masse,pk) -c=========================================================================== - -c ----------------------------------------------------------------- -c Inversion Nord/sud de la grille + interpollation sur les niveaux -c verticaux du modele. -c ----------------------------------------------------------------- - - implicit none - -#include "dimensions.h" -#include "paramet.h" -#include "comgeom2.h" -#include "comvert.h" -#include "comconst.h" -#include "guide.h" - - integer nlevnc - real psi(iip1,jjp1) - real u(iip1,jjp1,llm),v(iip1,jjm,llm) - real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm) - - real plnc(iip1,jjp1,nlevnc) - real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc) - real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1) - real qnc(iip1,jjp1,nlevnc) - - real zu(iip1,jjp1,llm),zv(iip1,jjm,llm) - real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm) - - real pext(iip1,jjp1,llm) - real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm) - real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm) - real plsnc(iip1,jjp1,llm) - - REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) - real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1) - real pkf(iip1,jjp1,llm) - real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm) - real prefkap,unskap - - - integer i,j,l - - -c ----------------------------------------------------------------- -c calcul de la pression au milieu des couches. -c ----------------------------------------------------------------- - - CALL pression( ip1jmp1, ap, bp, psi, p ) - call massdair(p,masse) - CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf) - -c .... Calcul de pls , pression au milieu des couches ,en Pascals - unskap=1./kappa - prefkap = preff ** kappa -c PRINT *,' Pref kappa unskap ',preff,kappa,unskap - DO l = 1, llm - DO j=1,jjp1 - DO i =1, iip1 - pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap - ENDDO - ENDDO - ENDDO - - -c ----------------------------------------------------------------- -c calcul des pressions pour les grilles u et v -c ----------------------------------------------------------------- - - do l=1,llm - do j=1,jjp1 - do i=1,iip1 - pext(i,j,l)=pls(i,j,l)*aire(i,j) - enddo - enddo - enddo - call massbar(pext, pbarx, pbary ) - do l=1,llm - do j=1,jjp1 - do i=1,iip1 - plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu(i,j) - plsnc(i,jjp1+1-j,l)=pls(i,j,l) - enddo - enddo - enddo - do l=1,llm - do j=1,jjm - do i=1,iip1 - plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev(i,j) - enddo - enddo - enddo - -c ----------------------------------------------------------------- - - if (guide_P) then - do j=1,jjp1 - do i=1,iim - ps(i,j)=psnc(i,jjp1+1-j) - enddo - ps(iip1,j)=ps(1,j) - enddo - endif - - -c ----------------------------------------------------------------- - call pres2lev(unc,zu,nlevnc,llm,plnc,plunc,iip1,jjp1,ok_invertp) - call pres2lev(vnc,zv,nlevnc,llm,plnc,plvnc,iip1,jjm,ok_invertp ) - call pres2lev(tnc,zt,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp) - call pres2lev(qnc,zq,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp) - -c call dump2d(iip1,jjp1,ps,'PS ') -c call dump2d(iip1,jjp1,psu,'PS ') -c call dump2d(iip1,jjm,psv,'PS ') -c Inversion Nord/Sud - do l=1,llm - do j=1,jjp1 - do i=1,iim - u(i,j,l)=zu(i,jjp1+1-j,l) - t(i,j,l)=zt(i,jjp1+1-j,l) - q(i,j,l)=zq(i,jjp1+1-j,l) - enddo - u(iip1,j,l)=u(1,j,l) - t(iip1,j,l)=t(1,j,l) - q(iip1,j,l)=q(1,j,l) - enddo - enddo - - do l=1,llm - do j=1,jjm - do i=1,iim - v(i,j,l)=zv(i,jjm+1-j,l) - enddo - v(iip1,j,l)=v(1,j,l) - enddo - enddo - - return - end - -c=========================================================================== - subroutine nat2gcm(u,v,t,rh,pk,ucov,vcov,teta,q) -c=========================================================================== - - implicit none - -#include "dimensions.h" -#include "paramet.h" -#include "comgeom2.h" -#include "comconst.h" -#include "comvert.h" -#include "guide.h" - - real u(iip1,jjp1,llm),v(iip1,jjm,llm) - real t(iip1,jjp1,llm),pk(iip1,jjp1,llm),rh(iip1,jjp1,llm) - real ps(iip1,jjp1) - - real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm) - real teta(iip1,jjp1,llm),q(iip1,jjp1,llm) - - real pres(iip1,jjp1,llm),qsat(iip1,jjp1,llm) - - real unskap - - integer i,j,l - - - print*,'Entree dans nat2gcm' -c ucov(:,:,:)=0. -c do l=1,llm -c ucov(:,2:jjm,l)=u(:,2:jjm,l)*cu(:,2:jjm) -c enddo -c ucov(iip1,:,:)=ucov(1,:,:) - -c teta(:,:,:)=t(:,:,:)*cpp/pk(:,:,:) -c teta(iip1,:,:)=teta(1,:,:) - -c calcul de ucov et de la temperature potentielle - do l=1,llm - do j=1,jjp1 - do i=1,iim - ucov(i,j,l)=u(i,j,l)*cu(i,j) - teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l) - enddo - ucov(iip1,j,l)=ucov(1,j,l) - teta(iip1,j,l)=teta(1,j,l) - enddo - do i=1,iip1 - ucov(i,1,l)=0. - ucov(i,jjp1,l)=0. - teta(i,1,l)=teta(1,1,l) - teta(i,jjp1,l)=teta(1,jjp1,l) - enddo - enddo - -c calcul de ucov - do l=1,llm - do j=1,jjm - do i=1,iim - vcov(i,j,l)=v(i,j,l)*cv(i,j) - enddo - vcov(iip1,j,l)=vcov(1,j,l) - enddo - enddo - -c call dump2d(iip1,jjp1,teta,'TETA EN BAS ') -c call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT ') - -c Humidite relative -> specifique -c ------------------------------- - if (guide_hr) then -c FINALEMENT ON GUIDE EN HUMIDITE RELATIVE - print*,'calcul de unskap' - unskap = 1./ kappa - print*,'calcul de pres' - pres(:,:,:)=preff*(pk(:,:,:)/cpp)**unskap - print*,'calcul de qsat' - call q_sat(iip1*jjp1*llm,t,pres,qsat) - print*,'calcul de q' -c ATTENTION : humidites relatives en % - rh(:,:,:)=max(rh(:,:,:)*0.01,1.e-6) - q(:,:,:)=qsat(:,:,:)*rh(:,:,:) - print*,'calcul de q OK' - call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1 ') - call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE ') - else - q(:,:,:)=rh(:,:,:) - print*,'calcul de q OK' - call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1 ') - endif - - return - end - - - -c=========================================================================== - subroutine correctbid(iim,nl,x) -c=========================================================================== - integer iim,nl - real x(iim+1,nl) - integer i,l - real zz - - do l=1,nl - do i=2,iim-1 - if(abs(x(i,l)).gt.1.e10) then - zz=0.5*(x(i-1,l)+x(i+1,l)) -c print*,'correction ',i,l,x(i,l),zz - x(i,l)=zz - endif - enddo - enddo - return - end