! ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/guide.F,v 1.3.4.1 2006/11/06 15:51:16 fairhead Exp $ ! MODULE guide_p_mod !======================================================================= ! Auteur: F.Hourdin ! F. Codron 01/09 !======================================================================= USE getparam USE Write_Field_p use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close IMPLICIT NONE ! --------------------------------------------- ! Declarations des cles logiques et parametres ! --------------------------------------------- INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav INTEGER, PRIVATE, SAVE :: nlevnc LOGICAL, PRIVATE, SAVE :: guide_u,guide_v,guide_T,guide_Q,guide_P LOGICAL, PRIVATE, SAVE :: guide_hr,guide_teta LOGICAL, PRIVATE, SAVE :: guide_BL,guide_reg,guide_add,gamma4,guide_zon LOGICAL, PRIVATE, SAVE :: guide_modele,invert_p,invert_y,ini_anal LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav REAL, PRIVATE, SAVE :: tau_min_u,tau_max_u REAL, PRIVATE, SAVE :: tau_min_v,tau_max_v REAL, PRIVATE, SAVE :: tau_min_T,tau_max_T REAL, PRIVATE, SAVE :: tau_min_Q,tau_max_Q REAL, PRIVATE, SAVE :: tau_min_P,tau_max_P REAL, PRIVATE, SAVE :: lat_min_g,lat_max_g REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g REAL, PRIVATE, SAVE :: tau_lon,tau_lat REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_T,alpha_Q REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P,alpha_pcor ! --------------------------------------------- ! Variables de guidage ! --------------------------------------------- ! Variables des fichiers de guidage REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: unat1,unat2 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: vnat1,vnat2 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: tnat1,tnat2 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: qnat1,qnat2 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: psnat1,psnat2 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc,bpnc ! Variables aux dimensions du modele REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: ugui1,ugui2 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: vgui1,vgui2 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tgui1,tgui2 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qgui1,qgui2 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1,psgui2 INTEGER,SAVE,PRIVATE :: ijb_u,ijb_v,ije_u,ije_v,ijn_u,ijn_v INTEGER,SAVE,PRIVATE :: jjb_u,jjb_v,jje_u,jje_v,jjn_u,jjn_v CONTAINS !======================================================================= SUBROUTINE guide_init IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "netcdf.inc" INCLUDE "control.h" INTEGER :: error,ncidpl,rid,rcod CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: modname = 'guide_init' ! --------------------------------------------- ! Lecture des parametres: ! --------------------------------------------- ! Variables guidees 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_teta',.false.,guide_teta,'guidage de T par Teta') CALL getpar('guide_add',.false.,guide_add,'for�age constant?') CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale') ! 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') CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim') ! Sauvegarde du for�age CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage') CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage') ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois. IF (iguide_sav.GT.0) THEN iguide_sav=day_step/iguide_sav ELSE iguide_sav=day_step*iguide_sav ENDIF ! Guidage regional seulement (sinon constant ou suivant le zoom) CALL getpar('guide_reg',.false.,guide_reg,'guidage regional') CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ') CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ') CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ') CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ') CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ') CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ') ! Parametres pour lecture des fichiers CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage') CALL getpar('iguide_int',4,iguide_int,'freq. lecture guidage') IF (iguide_int.GT.0) THEN iguide_int=day_step/iguide_int ELSE iguide_int=day_step*iguide_int ENDIF CALL getpar('guide_modele',.false.,guide_modele,'guidage niveaux modele') CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse') CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses') CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S') CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P') ! --------------------------------------------- ! Determination du nombre de niveaux verticaux ! des fichiers guidage ! --------------------------------------------- ncidpl=-99 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) elseif (guide_v) then if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl) elseif (guide_T) then if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl) elseif (guide_Q) then if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl) endif endif error=NF_INQ_DIMID(ncidpl,'LEVEL',rid) IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid) IF (error.NE.NF_NOERR) THEN print *,'Guide: probleme lecture niveaux pression' CALL abort_gcm(modname,abort_message,1) ENDIF error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc) print *,'Guide: nombre niveaux vert. nlevnc', nlevnc rcod = nf90_close(ncidpl) ! --------------------------------------------- ! Allocation des variables ! --------------------------------------------- abort_message='pb in allocation guide' ALLOCATE(apnc(nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(bpnc(nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) apnc=0.;bpnc=0. ALLOCATE(alpha_pcor(llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_u(ip1jmp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_v(ip1jm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_T(ip1jmp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_Q(ip1jmp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_P(ip1jmp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0 IF (guide_u) THEN ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(ugui1(ip1jmp1,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(ugui2(ip1jmp1,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) unat1=0.;unat2=0.;ugui1=0.;ugui2=0. ENDIF IF (guide_T) THEN ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(tgui1(ip1jmp1,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(tgui2(ip1jmp1,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0. ENDIF IF (guide_Q) THEN ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(qgui1(ip1jmp1,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(qgui2(ip1jmp1,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0. ENDIF IF (guide_v) THEN ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(vgui1(ip1jm,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(vgui2(ip1jm,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0. ENDIF IF (guide_P.OR.guide_modele) THEN ALLOCATE(psnat1(iip1,jjp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(psnat2(iip1,jjp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) psnat1=0.;psnat2=0.; ENDIF IF (guide_P) THEN ALLOCATE(psgui2(ip1jmp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(psgui1(ip1jmp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) psgui1=0.;psgui2=0. ENDIF ! --------------------------------------------- ! Lecture du premier etat de guidage. ! --------------------------------------------- IF (guide_2D) THEN CALL guide_read2D(1) ELSE CALL guide_read(1) ENDIF IF (guide_v) vnat1=vnat2 IF (guide_u) unat1=unat2 IF (guide_T) tnat1=tnat2 IF (guide_Q) qnat1=qnat2 IF (guide_P.OR.guide_modele) psnat1=psnat2 END SUBROUTINE guide_init !======================================================================= SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps) use parallel IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "control.h" INCLUDE "comconst.h" INCLUDE "comvert.h" ! Variables entree INTEGER, INTENT(IN) :: itau !pas de temps REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse REAL, DIMENSION (ip1jm,llm), INTENT(INOUT) :: vcov REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps ! Variables locales LOGICAL, SAVE :: first=.TRUE. LOGICAL :: f_out ! sortie guidage REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage REAL, DIMENSION (ip1jmp1,llm) :: p ! besoin si guide_P ! Compteurs temps: INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage REAL :: ditau, dday_step REAL :: tau,reste ! position entre 2 etats de guidage REAL, SAVE :: factt ! pas de temps en fraction de jour INTEGER :: l ijb_u=ij_begin ; ije_u=ij_end ; ijn_u=ije_u-ijb_u+1 jjb_u=jj_begin ; jje_u=jj_end ; jjn_u=jje_u-jjb_u+1 ijb_v=ij_begin ; ije_v=ij_end ; ijn_v=ije_v-ijb_v+1 jjb_v=jj_begin ; jje_v=jj_end ; jjn_v=jje_v-jjb_v+1 IF (pole_sud) THEN ije_v=ij_end-iip1 jje_v=jj_end-1 ijn_v=ije_v-ijb_v+1 jjn_v=jje_v-jjb_v+1 ENDIF PRINT *,'---> on rentre dans guide_main' ! CALL AllGather_Field(ucov,ip1jmp1,llm) ! CALL AllGather_Field(vcov,ip1jm,llm) ! CALL AllGather_Field(teta,ip1jmp1,llm) ! CALL AllGather_Field(ps,ip1jmp1,1) ! CALL AllGather_Field(q,ip1jmp1,llm) !----------------------------------------------------------------------- ! Initialisations au premier passage !----------------------------------------------------------------------- IF (first) THEN first=.FALSE. CALL guide_init itau_test=1001 step_rea=1 count_no_rea=0 ! Calcul des constantes de rappel factt=dtvr*iperiod/daysec 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) ! correction de rappel dans couche limite 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 ! ini_anal: etat initial egal au guidage IF (ini_anal) THEN CALL guide_interp(ps,teta) IF (guide_u) ucov(ijb_u:ije_u,:)=ugui2(ijb_u:ije_u,:) IF (guide_v) vcov(ijb_v:ije_v,:)=ugui2(ijb_v:ije_v,:) IF (guide_T) teta(ijb_u:ije_u,:)=tgui2(ijb_u:ije_u,:) IF (guide_Q) q(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:) IF (guide_P) THEN ps(ijb_u:ije_u)=psgui2(ijb_u:ije_u) CALL pression_p(ip1jmp1,ap,bp,ps,p) CALL massdair_p(p,masse) ENDIF RETURN ENDIF ! Verification structure guidage IF (guide_u) THEN CALL writefield_p('unat',unat1) CALL writefield_p('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/))) ENDIF IF (guide_T) THEN CALL writefield_p('tnat',tnat1) CALL writefield_p('teta',RESHAPE(teta,(/iip1,jjp1,llm/))) ENDIF ENDIF !first !----------------------------------------------------------------------- ! Lecture des fichiers de guidage ? !----------------------------------------------------------------------- IF (iguide_read.NE.0) THEN ditau=real(itau) dday_step=real(day_step) IF (iguide_read.LT.0) THEN tau=ditau/dday_step/FLOAT(iguide_read) ELSE tau=FLOAT(iguide_read)*ditau/dday_step ENDIF reste=tau-AINT(tau) IF (reste.EQ.0.) THEN IF (itau_test.EQ.itau) THEN write(*,*)'deuxieme passage de advreel a itau=',itau stop ELSE IF (guide_v) vnat1(:,jjb_v:jje_v,:)=vnat2(:,jjb_v:jje_v,:) IF (guide_u) unat1(:,jjb_u:jje_u,:)=unat2(:,jjb_u:jje_u,:) IF (guide_T) tnat1(:,jjb_u:jje_u,:)=tnat2(:,jjb_u:jje_u,:) IF (guide_Q) qnat1(:,jjb_u:jje_u,:)=qnat2(:,jjb_u:jje_u,:) IF (guide_P.OR.guide_modele) psnat1(:,jjb_u:jje_u)=psnat2(:,jjb_u:jje_u) step_rea=step_rea+1 itau_test=itau print*,'Lecture fichiers guidage, pas ',step_rea, & 'apres ',count_no_rea,' non lectures' IF (guide_2D) THEN CALL guide_read2D(step_rea) ELSE CALL guide_read(step_rea) ENDIF count_no_rea=0 ENDIF ELSE count_no_rea=count_no_rea+1 ENDIF ENDIF !iguide_read=0 !----------------------------------------------------------------------- ! Interpolation et conversion des champs de guidage !----------------------------------------------------------------------- IF (MOD(itau,iguide_int).EQ.0) THEN CALL guide_interp(ps,teta) ENDIF ! Repartition entre 2 etats de guidage IF (iguide_read.NE.0) THEN tau=reste ELSE tau=1. ENDIF !----------------------------------------------------------------------- ! Ajout des champs de guidage !----------------------------------------------------------------------- ! Sauvegarde du guidage? f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav) IF (f_out) CALL guide_out("S",jjp1,1,ps) if (guide_u) then if (guide_add) then f_add(ijb_u:ije_u,:)=(1.-tau)*ugui1(ijb_u:ije_u,:)+tau*ugui2(ijb_u:ije_u,:) else f_add(ijb_u:ije_u,:)=(1.-tau)*ugui1(ijb_u:ije_u,:)+tau*ugui2(ijb_u:ije_u,:)-ucov(ijb_u:ije_u,:) endif if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add) CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u) IF (f_out) CALL guide_out("U",jjp1,llm,f_add/factt) ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:) endif if (guide_T) then if (guide_add) then f_add(ijb_u:ije_u,:)=(1.-tau)*tgui1(ijb_u:ije_u,:)+tau*tgui2(ijb_u:ije_u,:) else f_add(ijb_u:ije_u,:)=(1.-tau)*tgui1(ijb_u:ije_u,:)+tau*tgui2(ijb_u:ije_u,:)-teta(ijb_u:ije_u,:) endif if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add) CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T) IF (f_out) CALL guide_out("T",jjp1,llm,f_add/factt) teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:) endif if (guide_P) then if (guide_add) then f_add(ijb_u:ije_u,1)=(1.-tau)*psgui1(ijb_u:ije_u)+tau*psgui2(ijb_u:ije_u) else f_add(ijb_u:ije_u,1)=(1.-tau)*psgui1(ijb_u:ije_u)+tau*psgui2(ijb_u:ije_u)-ps(ijb_u:ije_u) endif if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1)) CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P) IF (f_out) CALL guide_out("P",jjp1,1,f_add(1:ip1jmp1,1)/factt) ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1) CALL pression_p(ip1jmp1,ap,bp,ps,p) CALL massdair_p(p,masse) endif if (guide_Q) then if (guide_add) then f_add(ijb_u:ije_u,:)=(1.-tau)*qgui1(ijb_u:ije_u,:)+tau*qgui2(ijb_u:ije_u,:) else f_add(ijb_u:ije_u,:)=(1.-tau)*qgui1(ijb_u:ije_u,:)+tau*qgui2(ijb_u:ije_u,:)-q(ijb_u:ije_u,:) endif if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add) CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q) IF (f_out) CALL guide_out("Q",jjp1,llm,f_add/factt) q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:) endif if (guide_v) then if (guide_add) then f_add(ijb_v:ije_v,:)=(1.-tau)*vgui1(ijb_v:ije_v,:)+tau*vgui2(ijb_v:ije_v,:) else f_add(ijb_v:ije_v,:)=(1.-tau)*vgui1(ijb_v:ije_v,:)+tau*vgui2(ijb_v:ije_v,:)-vcov(ijb_v:ije_v,:) endif if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:)) CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v) IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:)/factt) vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:) endif END SUBROUTINE guide_main !======================================================================= SUBROUTINE guide_addfield(hsize,vsize,field,alpha) ! field1=a*field1+alpha*field2 IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" ! input variables INTEGER, INTENT(IN) :: hsize INTEGER, INTENT(IN) :: vsize REAL, DIMENSION(hsize), INTENT(IN) :: alpha REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field ! Local variables INTEGER :: l IF (hsize==ip1jm) THEN do l=1,vsize field(ijb_v:ije_v,l)=alpha(ijb_v:ije_v)*field(ijb_v:ije_v,l)*alpha_pcor(l) enddo ELSE do l=1,vsize field(ijb_u:ije_u,l)=alpha(ijb_u:ije_u)*field(ijb_u:ije_u,l)*alpha_pcor(l) enddo ENDIF END SUBROUTINE guide_addfield !======================================================================= SUBROUTINE guide_zonave(typ,hsize,vsize,field) IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "comgeom.h" INCLUDE "comconst.h" ! input/output variables INTEGER, INTENT(IN) :: typ INTEGER, INTENT(IN) :: vsize INTEGER, INTENT(IN) :: hsize REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field ! Local variables LOGICAL, SAVE :: first=.TRUE. INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain INTEGER :: i,j,l,ij REAL, DIMENSION (iip1) :: lond ! longitude in Deg. REAL, DIMENSION (hsize,vsize):: fieldm ! zon-averaged field IF (first) THEN first=.FALSE. !Compute domain for averaging lond=rlonu*180./pi imin(1)=1;imax(1)=iip1; imin(2)=1;imax(2)=iip1; IF (guide_reg) THEN DO i=1,iim IF (lond(i).LT.lon_min_g) imin(1)=i IF (lond(i).LE.lon_max_g) imax(1)=i ENDDO lond=rlonv*180./pi DO i=1,iim IF (lond(i).LT.lon_min_g) imin(2)=i IF (lond(i).LE.lon_max_g) imax(2)=i ENDDO ENDIF ENDIF fieldm=0. IF (hsize==jjm) THEN DO l=1,vsize ! Compute zonal average DO j=jjb_v,jje_v DO i=imin(typ),imax(typ) ij=(j-1)*iip1+i fieldm(j,l)=fieldm(j,l)+field(ij,l) ENDDO ENDDO fieldm(:,l)=fieldm(:,l)/FLOAT(imax(typ)-imin(typ)+1) ! Compute forcing DO j=jjb_v,jje_v DO i=1,iip1 ij=(j-1)*iip1+i field(ij,l)=fieldm(j,l) ENDDO ENDDO ENDDO ELSE DO l=1,vsize ! Compute zonal average DO j=jjb_v,jje_v DO i=imin(typ),imax(typ) ij=(j-1)*iip1+i fieldm(j,l)=fieldm(j,l)+field(ij,l) ENDDO ENDDO fieldm(:,l)=fieldm(:,l)/FLOAT(imax(typ)-imin(typ)+1) ! Compute forcing DO j=jjb_u,jje_u DO i=1,iip1 ij=(j-1)*iip1+i field(ij,l)=fieldm(j,l) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE guide_zonave !======================================================================= SUBROUTINE guide_interp(psi,teta) USE parallel USE mod_hallo USE Bands IMPLICIT NONE include "dimensions.h" include "paramet.h" include "comvert.h" include "comgeom2.h" include "comconst.h" REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm LOGICAL, SAVE :: first=.TRUE. ! Variables pour niveaux pression: REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire REAL, DIMENSION (iip1,jjp1,llm) :: pbarx REAL, DIMENSION (iip1,jjm,llm) :: pbary ! Variables pour fonction Exner (P milieu couche) REAL, DIMENSION (iip1,jjp1,llm) :: pk, pkf REAL, DIMENSION (iip1,jjp1,llm) :: alpha, beta REAL, DIMENSION (iip1,jjp1) :: pks REAL :: prefkap,unskap ! Pression de vapeur saturante REAL, DIMENSION (ip1jmp1,llm) :: qsat !Variables intermediaires interpolation REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2 REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2 INTEGER :: i,j,l,ij TYPE(Request) :: Req print *,'Guide: conversion variables guidage' ! ----------------------------------------------------------------- ! Calcul des niveaux de pression champs guidage ! ----------------------------------------------------------------- if (guide_modele) then do i=1,iip1 do j=jjb_u,jje_u do l=1,nlevnc plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j) plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j) enddo enddo enddo else do i=1,iip1 do j=jjb_u,jje_u do l=1,nlevnc plnc2(i,j,l)=apnc(l) plnc1(i,j,l)=apnc(l) enddo enddo enddo endif if (first) then first=.FALSE. print*,'Guide: verification ordre niveaux verticaux' print*,'LMDZ :' do l=1,llm print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. & +psi(1,jje_u)*(bp(l)+bp(l+1))/2. enddo print*,'Fichiers guidage' do l=1,nlevnc print*,'PL(',l,')=',plnc2(1,jjb_u,l) enddo print *,'inversion de l''ordre: invert_p=',invert_p if (guide_u) then do l=1,nlevnc print*,'U(',l,')=',unat2(1,jjb_u,l) enddo endif if (guide_T) then do l=1,nlevnc print*,'T(',l,')=',tnat2(1,jjb_u,l) enddo endif endif ! ----------------------------------------------------------------- ! Calcul niveaux pression modele ! ----------------------------------------------------------------- CALL pression_p( ip1jmp1, ap, bp, psi, p ) CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf) ! .... Calcul de pls , pression au milieu des couches ,en Pascals unskap=1./kappa prefkap = preff ** kappa DO l = 1, llm DO j=jjb_u,jje_u DO i =1, iip1 pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap ENDDO ENDDO ENDDO ! calcul des pressions pour les grilles u et v do l=1,llm do j=jjb_u,jje_u do i=1,iip1 pext(i,j,l)=pls(i,j,l)*aire(i,j) enddo enddo enddo CALL Register_SwapFieldHallo(pext,pext,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req) CALL SendRequest(Req) CALL WaitRequest(Req) call massbar_p(pext, pbarx, pbary ) do l=1,llm do j=jjb_u,jje_u do i=1,iip1 plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j) plsnc(i,j,l)=pls(i,j,l) enddo enddo enddo do l=1,llm do j=jjb_v,jje_v do i=1,iip1 plvnc(i,j,l)=pbary(i,j,l)/airev(i,j) enddo enddo enddo ! ----------------------------------------------------------------- ! Interpolation champs guidage sur niveaux modele (+inversion N/S) ! Conversion en variables gcm (ucov, vcov...) ! ----------------------------------------------------------------- if (guide_P) then do j=jjb_u,jje_u do i=1,iim ij=(j-1)*iip1+i psgui1(ij)=psnat1(i,j) psgui2(ij)=psnat2(i,j) enddo psgui1(iip1*j)=psnat1(1,j) psgui2(iip1*j)=psnat2(1,j) enddo endif IF (guide_u) THEN CALL pres2lev(unat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, & plnc1(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p) CALL pres2lev(unat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, & plnc2(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p) do l=1,llm do j=jjb_u,jje_u do i=1,iim ij=(j-1)*iip1+i ugui1(ij,l)=zu1(i,j,l)*cu(i,j) ugui2(ij,l)=zu2(i,j,l)*cu(i,j) enddo ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l) ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l) enddo do i=1,iip1 ugui1(i,l)=0. ugui1(ip1jm+i,l)=0. ugui2(i,l)=0. ugui2(ip1jm+i,l)=0. enddo enddo ENDIF IF (guide_T) THEN CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, & plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p) CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, & plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p) do l=1,llm do j=jjb_u,jje_u IF (guide_teta) THEN do i=1,iim ij=(j-1)*iip1+i tgui1(ij,l)=zu1(i,j,l) tgui2(ij,l)=zu2(i,j,l) enddo ELSE do i=1,iim ij=(j-1)*iip1+i tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l) tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l) enddo ENDIF tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l) tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l) enddo do i=1,iip1 tgui1(i,l)=tgui1(1,l) tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l) tgui2(i,l)=tgui2(1,l) tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l) enddo enddo ENDIF IF (guide_v) THEN CALL pres2lev(vnat1(:,jjb_v:jje_v,:),zv1(:,jjb_v:jje_v,:),nlevnc,llm, & plnc1(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p) CALL pres2lev(vnat2(:,jjb_v:jje_v,:),zv2(:,jjb_v:jje_v,:),nlevnc,llm, & plnc2(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p) do l=1,llm do j=jjb_v,jje_v do i=1,iim ij=(j-1)*iip1+i vgui1(ij,l)=zv1(i,j,l)*cv(i,j) vgui2(ij,l)=zv2(i,j,l)*cv(i,j) enddo vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l) vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l) enddo enddo ENDIF IF (guide_Q) THEN ! On suppose qu'on a la bonne variable dans le fichier de guidage: ! Hum.Rel si guide_hr, Hum.Spec. sinon. CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, & plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p) CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, & plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p) do l=1,llm do j=jjb_u,jjb_v do i=1,iim ij=(j-1)*iip1+i qgui1(ij,l)=zu1(i,j,l) qgui2(ij,l)=zu2(i,j,l) enddo qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l) qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l) enddo do i=1,iip1 qgui1(i,l)=qgui1(1,l) qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l) qgui2(i,l)=qgui2(1,l) qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l) enddo enddo IF (guide_hr) THEN CALL q_sat(iip1*jjn_u*llm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp, & plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:)) qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01 !hum. rel. en % qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01 ENDIF ENDIF END SUBROUTINE guide_interp !======================================================================= SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha) ! Calcul des constantes de rappel alpha (=1/tau) implicit none include "dimensions.h" include "paramet.h" include "comconst.h" include "comgeom2.h" include "serre.h" ! input arguments : INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1) INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour REAL, INTENT(IN) :: taumin,taumax ! output arguments: REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha ! local variables: LOGICAL, SAVE :: first=.TRUE. REAL, SAVE :: gamma,dxdy_min,dxdy_max REAL, DIMENSION (iip1,jjp1) :: zdx,zdy REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu REAL, DIMENSION (iip1,jjm) :: dxdyv real dxdy_ real zlat,zlon real alphamin,alphamax,xi integer i,j,ilon,ilat alphamin=factt/taumax alphamax=factt/taumin IF (guide_reg.OR.guide_add) THEN alpha=alphamax !----------------------------------------------------------------------- ! guide_reg: alpha=alpha_min dans region, 0. sinon. !----------------------------------------------------------------------- IF (guide_reg) THEN do j=1,pjm do i=1,pim if (typ.eq.2) then zlat=rlatu(j)*180./pi zlon=rlonu(i)*180./pi elseif (typ.eq.1) then zlat=rlatu(j)*180./pi zlon=rlonv(i)*180./pi elseif (typ.eq.3) then zlat=rlatv(j)*180./pi zlon=rlonv(i)*180./pi endif alpha(i,j)=alphamax/16.* & (1.+tanh((zlat-lat_min_g)/tau_lat))* & (1.+tanh((lat_max_g-zlat)/tau_lat))* & (1.+tanh((zlon-lon_min_g)/tau_lon))* & (1.+tanh((lon_max_g-zlon)/tau_lon)) enddo enddo ENDIF ELSE !----------------------------------------------------------------------- ! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom. !----------------------------------------------------------------------- !Calcul de l'aire des mailles 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 IF (typ.EQ.2) THEN 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 ENDIF IF (typ.EQ.3) THEN do j=1,jjm do i=1,iip1 dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1)) enddo enddo ENDIF ! Premier appel: calcul des aires min et max et de gamma. IF (first) THEN first=.FALSE. ! coordonnees du centre du zoom CALL coordij(clon,clat,ilon,ilat) ! aire de la maille au centre du zoom dxdy_min=dxdys(ilon,ilat) ! 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 ! Calcul de gamma 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 !first do j=1,pjm do i=1,pim if (typ.eq.1) then dxdy_=dxdys(i,j) zlat=rlatu(j)*180./pi elseif (typ.eq.2) then dxdy_=dxdyu(i,j) zlat=rlatu(j)*180./pi elseif (typ.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 ! 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_g.le.zlat .and. zlat.le.lat_max_g) then alpha(i,j)=xi*alphamin+(1.-xi)*alphamax else alpha(i,j)=0. endif endif enddo enddo ENDIF ! guide_reg END SUBROUTINE tau2alpha !======================================================================= SUBROUTINE guide_read(timestep) IMPLICIT NONE #include "netcdf.inc" #include "dimensions.h" #include "paramet.h" INTEGER, INTENT(IN) :: timestep LOGICAL, SAVE :: first=.TRUE. ! Identification fichiers et variables NetCDF: INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps INTEGER :: ncidpl,varidpl,varidap,varidbp ! Variables auxiliaires NetCDF: INTEGER, DIMENSION(4) :: start,count INTEGER :: status,rcode ! ----------------------------------------------------------------- ! Premier appel: initialisation de la lecture des fichiers ! ----------------------------------------------------------------- if (first) then ncidpl=-99 print*,'Guide: ouverture des fichiers guidage ' ! Niveaux de pression si non constants if (guide_modele) then print *,'Lecture du guidage 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 ! 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 ! 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 ! 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 ! 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 ! 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 ! Coordonnee verticale if (.not.guide_modele) then rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl) IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl) print*,'ncidpl,varidpl',ncidpl,varidpl endif ! Coefs ap, bp pour calcul de la pression aux differents niveaux 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 apnc=apnc*100.! conversion en Pascals bpnc(:)=0. endif first=.FALSE. endif ! (first) ! ----------------------------------------------------------------- ! lecture des champs u, v, T, Q, ps ! ----------------------------------------------------------------- ! dimensions pour les champs scalaires et le vent zonal start(1)=1 start(2)=1 start(3)=1 start(4)=timestep count(1)=iip1 count(2)=jjp1 count(3)=nlevnc count(4)=1 ! Vent zonal if (guide_u) then #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2) #else status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2) #endif IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,unat2) ENDIF endif ! Temperature if (guide_T) then #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2) #else status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2) #endif IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,tnat2) ENDIF endif ! Humidite if (guide_Q) then #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2) #else status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2) #endif IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,qnat2) ENDIF endif ! Vent meridien if (guide_v) then count(2)=jjm #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2) #else status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2) #endif IF (invert_y) THEN CALL invert_lat(iip1,jjm,nlevnc,vnat2) ENDIF endif ! Pression de surface if ((guide_P).OR.(guide_modele)) then start(3)=timestep start(4)=0 count(2)=jjp1 count(3)=1 count(4)=0 #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2) #else status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2) #endif IF (invert_y) THEN CALL invert_lat(iip1,jjp1,1,psnat2) ENDIF endif END SUBROUTINE guide_read !======================================================================= SUBROUTINE guide_read2D(timestep) IMPLICIT NONE #include "netcdf.inc" #include "dimensions.h" #include "paramet.h" INTEGER, INTENT(IN) :: timestep LOGICAL, SAVE :: first=.TRUE. ! Identification fichiers et variables NetCDF: INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps INTEGER :: ncidpl,varidpl,varidap,varidbp ! Variables auxiliaires NetCDF: INTEGER, DIMENSION(4) :: start,count INTEGER :: status,rcode ! Variables for 3D extension: REAL, DIMENSION (jjp1,llm) :: zu REAL, DIMENSION (jjm,llm) :: zv INTEGER :: i ! ----------------------------------------------------------------- ! Premier appel: initialisation de la lecture des fichiers ! ----------------------------------------------------------------- if (first) then ncidpl=-99 print*,'Guide: ouverture des fichiers guidage ' ! Niveaux de pression si non constants if (guide_modele) then print *,'Lecture du guidage 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 ! 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 ! 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 ! 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 ! 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 ! 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 ! Coordonnee verticale if (.not.guide_modele) then rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl) IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl) print*,'ncidpl,varidpl',ncidpl,varidpl endif ! Coefs ap, bp pour calcul de la pression aux differents niveaux 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 apnc=apnc*100.! conversion en Pascals bpnc(:)=0. endif first=.FALSE. endif ! (first) ! ----------------------------------------------------------------- ! lecture des champs u, v, T, Q, ps ! ----------------------------------------------------------------- ! dimensions pour les champs scalaires et le vent zonal start(1)=1 start(2)=1 start(3)=1 start(4)=timestep count(1)=1 count(2)=jjp1 count(3)=nlevnc count(4)=1 ! Vent zonal if (guide_u) then #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu) #else status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu) #endif DO i=1,iip1 unat2(i,:,:)=zu(:,:) ENDDO IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,unat2) ENDIF endif ! Temperature if (guide_T) then #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu) #else status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu) #endif DO i=1,iip1 tnat2(i,:,:)=zu(:,:) ENDDO IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,tnat2) ENDIF endif ! Humidite if (guide_Q) then #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu) #else status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu) #endif DO i=1,iip1 qnat2(i,:,:)=zu(:,:) ENDDO IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,qnat2) ENDIF endif ! Vent meridien if (guide_v) then count(2)=jjm #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv) #else status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv) #endif DO i=1,iip1 vnat2(i,:,:)=zv(:,:) ENDDO IF (invert_y) THEN CALL invert_lat(iip1,jjm,nlevnc,vnat2) ENDIF endif ! Pression de surface if ((guide_P).OR.(guide_modele)) then start(3)=timestep start(4)=0 count(2)=jjp1 count(3)=1 count(4)=0 #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1)) #else status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1)) #endif DO i=1,iip1 psnat2(i,:)=zu(:,1) ENDDO IF (invert_y) THEN CALL invert_lat(iip1,jjp1,1,psnat2) ENDIF endif END SUBROUTINE guide_read2D !======================================================================= SUBROUTINE guide_out(varname,hsize,vsize,field) USE parallel IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "netcdf.inc" INCLUDE "comgeom2.h" INCLUDE "comconst.h" INCLUDE "comvert.h" ! Variables entree CHARACTER, INTENT(IN) :: varname INTEGER, INTENT (IN) :: hsize,vsize REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field ! Variables locales INTEGER, SAVE :: timestep=0 ! Identites fichier netcdf INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev INTEGER, DIMENSION (3) :: dim3 INTEGER, DIMENSION (4) :: dim4,count,start INTEGER :: ierr, varid CALL gather_field(field,iip1*hsize,vsize,0) IF (mpi_rank /= 0) RETURN print *,'Guide: output timestep',timestep,'var ',varname IF (timestep.EQ.0) THEN ! ---------------------------------------------- ! initialisation fichier de sortie ! ---------------------------------------------- ! Ouverture du fichier ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid) ! Definition des dimensions ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu) ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv) ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu) ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv) ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev) ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim) ! Creation des variables dimensions ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu) ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv) ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu) ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv) ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev) ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu) ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv) ierr=NF_ENDDEF(nid) ! Enregistrement des variables dimensions #ifdef NC_DOUBLE ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi) ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi) ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi) ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi) ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs) ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu) ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv) #else ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi) ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi) ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi) ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi) ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs) ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu) ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv) #endif ! -------------------------------------------------------------------- ! Cr�ation des variables sauvegard�es ! -------------------------------------------------------------------- ierr = NF_REDEF(nid) ! Surface pressure (GCM) dim3=(/id_lonv,id_latu,id_tim/) ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid) ! Surface pressure (guidage) IF (guide_P) THEN dim3=(/id_lonv,id_latu,id_tim/) ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid) ENDIF ! Zonal wind IF (guide_u) THEN dim4=(/id_lonu,id_latu,id_lev,id_tim/) ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid) ENDIF ! Merid. wind IF (guide_v) THEN dim4=(/id_lonv,id_latv,id_lev,id_tim/) ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid) ENDIF ! Pot. Temperature IF (guide_T) THEN dim4=(/id_lonv,id_latu,id_lev,id_tim/) ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid) ENDIF ! Specific Humidity IF (guide_Q) THEN dim4=(/id_lonv,id_latu,id_lev,id_tim/) ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid) ENDIF ierr = NF_ENDDEF(nid) ierr = NF_CLOSE(nid) ENDIF ! timestep=0 ! -------------------------------------------------------------------- ! Enregistrement du champ ! -------------------------------------------------------------------- ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid) SELECT CASE (varname) CASE ("S") timestep=timestep+1 ierr = NF_INQ_VARID(nid,"SP",varid) start=(/1,1,timestep,0/) count=(/iip1,jjp1,1,0/) #ifdef NC_DOUBLE ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field) #else ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field) #endif CASE ("P") ierr = NF_INQ_VARID(nid,"ps",varid) start=(/1,1,timestep,0/) count=(/iip1,jjp1,1,0/) #ifdef NC_DOUBLE ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field) #else ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field) #endif CASE ("U") ierr = NF_INQ_VARID(nid,"ucov",varid) start=(/1,1,1,timestep/) count=(/iip1,jjp1,llm,1/) #ifdef NC_DOUBLE ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field) #else ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field) #endif CASE ("V") ierr = NF_INQ_VARID(nid,"vcov",varid) start=(/1,1,1,timestep/) count=(/iip1,jjm,llm,1/) #ifdef NC_DOUBLE ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field) #else ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field) #endif CASE ("T") ierr = NF_INQ_VARID(nid,"teta",varid) start=(/1,1,1,timestep/) count=(/iip1,jjp1,llm,1/) #ifdef NC_DOUBLE ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field) #else ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field) #endif CASE ("Q") ierr = NF_INQ_VARID(nid,"q",varid) start=(/1,1,1,timestep/) count=(/iip1,jjp1,llm,1/) #ifdef NC_DOUBLE ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field) #else ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field) #endif END SELECT ierr = NF_CLOSE(nid) END SUBROUTINE guide_out !=========================================================================== subroutine correctbid(iim,nl,x) 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)) print*,'correction ',i,l,x(i,l),zz x(i,l)=zz endif enddo enddo return end subroutine correctbid !=========================================================================== END MODULE guide_p_mod