! ! $Id$ ! MODULE guide_mod !======================================================================= ! Auteur: F.Hourdin ! F. Codron 01/09 !======================================================================= USE iniprint_mod_h USE getparam, only: ini_getparam, fin_getparam, getpar USE Write_Field use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, & nf90_inq_dimid, nf90_inquire_dimension use pres2lev_mod, only: pres2lev USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, & nf90_inq_dimid, nf90_inquire_dimension, nf90_float, nf90_def_var, & nf90_create, nf90_def_dim, nf90_open, nf90_unlimited, nf90_write, nf90_enddef, nf90_redef, & nf90_close, nf90_inq_varid, nf90_get_var, nf90_noerr, nf90_clobber, & nf90_64bit_offset, nf90_inq_dimid, nf90_inquire_dimension, nf90_put_var USE paramet_mod_h IMPLICIT NONE ! --------------------------------------------- ! Declarations des cles logiques et parametres ! --------------------------------------------- INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav INTEGER, PRIVATE, SAVE :: nlevnc, guide_plevs 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 :: invert_p,invert_y,ini_anal LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav,guide_modele !FC LOGICAL, PRIVATE, SAVE :: convert_Pa 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, PRIVATE, SAVE :: plim_guide_BL 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 :: pnat1,pnat2 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 CONTAINS !======================================================================= SUBROUTINE guide_init use netcdf, only: nf90_noerr USE control_mod, ONLY: day_step USE serre_mod, ONLY: grossismx USE dimensions_mod, ONLY: iim, jjm, llm, ndm IMPLICIT NONE INTEGER :: error,ncidpl,rid,rcod CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: modname = 'guide_init' CHARACTER (len = 20) :: namedim ! --------------------------------------------- ! Lecture des parametres: ! --------------------------------------------- call ini_getparam("nudging_parameters_out.txt") ! 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') if (guide_zon .and. abs(grossismx - 1.) > 0.01) & call abort_gcm("guide_init", & "zonal nudging requires grid regular in longitude", 1) ! 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') CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value') ! 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 if (iguide_sav == 0) then iguide_sav = huge(0) 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. interpolation vert') IF (iguide_int.EQ.0) THEN iguide_int=1 ELSEIF (iguide_int.GT.0) THEN iguide_int=day_step/iguide_int ELSE iguide_int=day_step*iguide_int ENDIF CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage') ! Pour compatibilite avec ancienne version avec guide_modele CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol') IF (guide_modele) THEN guide_plevs=1 ENDIF !FC CALL getpar('convert_Pa',.true.,convert_Pa,'Convert Pressure levels in Pa') ! Fin raccord 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') call fin_getparam ! --------------------------------------------- ! Determination du nombre de niveaux verticaux ! des fichiers guidage ! --------------------------------------------- ncidpl=-99 if (guide_plevs.EQ.1) then if (ncidpl.eq.-99) then rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl) if (rcod.NE.NF90_NOERR) THEN abort_message=' Nudging error -> no file apbp.nc' CALL abort_gcm(modname,abort_message,1) endif endif elseif (guide_plevs.EQ.2) then if (ncidpl.EQ.-99) then rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl) if (rcod.NE.NF90_NOERR) THEN abort_message=' Nudging error -> no file P.nc' CALL abort_gcm(modname,abort_message,1) endif endif elseif (guide_u) then if (ncidpl.eq.-99) then rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl) if (rcod.NE.NF90_NOERR) THEN CALL abort_gcm(modname, & ' Nudging error -> no file u.nc',1) endif endif elseif (guide_v) then if (ncidpl.eq.-99) then rcod=nf90_open('v.nc',nf90_nowrite,ncidpl) if (rcod.NE.NF90_NOERR) THEN CALL abort_gcm(modname, & ' Nudging error -> no file v.nc',1) endif endif elseif (guide_T) then if (ncidpl.eq.-99) then rcod=nf90_open('T.nc',nf90_nowrite,ncidpl) if (rcod.NE.NF90_NOERR) THEN CALL abort_gcm(modname, & ' Nudging error -> no file T.nc',1) endif endif elseif (guide_Q) then if (ncidpl.eq.-99) then rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl) if (rcod.NE.NF90_NOERR) THEN CALL abort_gcm(modname, & ' Nudging error -> no file hur.nc',1) endif endif endif error=nf90_inq_dimid(ncidpl,'LEVEL',rid) IF (error.NE.NF90_NOERR) error=nf90_inq_dimid(ncidpl,'PRESSURE',rid) IF (error.NE.NF90_NOERR) THEN CALL abort_gcm(modname,'Nudging: error reading pressure levels',1) ENDIF error=nf90_inquire_dimension(ncidpl,rid,len=nlevnc) write(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc rcod = nf90_close(ncidpl) ! --------------------------------------------- ! Allocation des variables ! --------------------------------------------- abort_message='nudging allocation error' 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(iip1, jjp1), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_Q(iip1, jjp1), 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_plevs.EQ.2) THEN ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) pnat1=0.;pnat2=0.; ENDIF IF (guide_P.OR.guide_plevs.EQ.1) 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_plevs.EQ.2) pnat1=pnat2 IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2 END SUBROUTINE guide_init !======================================================================= SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps) USE exner_hyb_m, ONLY: exner_hyb USE exner_milieu_m, ONLY: exner_milieu USE control_mod, ONLY: day_step, iperiod USE comconst_mod, ONLY: cpp, dtvr, daysec,kappa USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h IMPLICIT NONE ! 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 :: pk(ip1jmp1,llm) ! Exner at mid-layers REAL :: pks(ip1jmp1) ! Exner at the surface REAL :: unskap ! 1./kappa REAL, DIMENSION (ip1jmp1,llmp1) :: p ! Pressure at inter-layers ! 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 CHARACTER(LEN=20) :: modname="guide_main" CHARACTER (len = 80) :: abort_message !----------------------------------------------------------------------- ! 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(((plim_guide_BL-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=ugui2 IF (guide_v) vcov=ugui2 IF (guide_T) teta=tgui2 IF (guide_Q) q=qgui2 IF (guide_P) THEN ps=psgui2 CALL pression(ip1jmp1,ap,bp,ps,p) CALL massdair(p,masse) ENDIF RETURN ENDIF ! Verification structure guidage ! IF (guide_u) THEN ! CALL writefield('unat',unat1) ! CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/))) ! ENDIF ! IF (guide_T) THEN ! CALL writefield('tnat',tnat1) ! CALL writefield('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/REAL(iguide_read) ELSE tau=REAL(iguide_read)*ditau/dday_step ENDIF reste=tau-AINT(tau) IF (reste.EQ.0.) THEN IF (itau_test.EQ.itau) THEN write(lunout,*)trim(modname)//' second pass in advreel at itau=',& itau abort_message='stopped' CALL abort_gcm(modname,abort_message,1) ELSE IF (guide_v) vnat1=vnat2 IF (guide_u) unat1=unat2 IF (guide_T) tnat1=tnat2 IF (guide_Q) qnat1=qnat2 IF (guide_plevs.EQ.2) pnat1=pnat2 IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2 step_rea=step_rea+1 itau_test=itau write(*,*)trim(modname)//' Reading nudging files, step ',& step_rea,'after ',count_no_rea,' skips' 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) THEN ! compute pressures at layer interfaces CALL pression(ip1jmp1,ap,bp,ps,p) if (pressure_exner) then call exner_hyb(ip1jmp1,ps,p,pks,pk) else call exner_milieu(ip1jmp1,ps,p,pks,pk) endif unskap=1./kappa ! Now compute pressures at mid-layer do l=1,llm p(:,l)=preff*(pk(:,l)/cpp)**unskap enddo CALL guide_out("SP",jjp1,llm,p(:,1:llm)) ENDIF if (guide_u) then if (guide_add) then f_add=(1.-tau)*ugui1+tau*ugui2 else f_add=(1.-tau)*ugui1+tau*ugui2-ucov 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("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2) IF (f_out) CALL guide_out("u",jjp1,llm,ucov) IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt) ucov=ucov+f_add endif if (guide_T) then if (guide_add) then f_add=(1.-tau)*tgui1+tau*tgui2 else f_add=(1.-tau)*tgui1+tau*tgui2-teta 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("teta",jjp1,llm,f_add/factt) teta=teta+f_add endif if (guide_P) then if (guide_add) then f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2 else f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2-ps 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("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt) ps=ps+f_add(1:ip1jmp1,1) CALL pression(ip1jmp1,ap,bp,ps,p) CALL massdair(p,masse) endif if (guide_Q) then if (guide_add) then f_add=(1.-tau)*qgui1+tau*qgui2 else f_add=(1.-tau)*qgui1+tau*qgui2-q 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=q+f_add endif if (guide_v) then if (guide_add) then f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2 else f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2-vcov 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,vcov) IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2) IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt) vcov=vcov+f_add(1:ip1jm,:) endif END SUBROUTINE guide_main !======================================================================= SUBROUTINE guide_addfield(hsize,vsize,field,alpha) ! field1=a*field1+alpha*field2 IMPLICIT NONE ! 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 do l=1,vsize field(:,l)=alpha*field(:,l)*alpha_pcor(l) enddo END SUBROUTINE guide_addfield !======================================================================= SUBROUTINE guide_zonave(typ,hsize,vsize,field) USE comconst_mod, ONLY: pi USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h USE comgeom_mod_h IMPLICIT NONE ! 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. DO l=1,vsize ! Compute zonal average DO j=1,hsize 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)/REAL(imax(typ)-imin(typ)+1) ! Compute forcing DO j=1,hsize DO i=1,iip1 ij=(j-1)*iip1+i field(ij,l)=fieldm(j,l) ENDDO ENDDO ENDDO END SUBROUTINE guide_zonave !======================================================================= SUBROUTINE guide_interp(psi,teta) use exner_hyb_m, only: exner_hyb use exner_milieu_m, only: exner_milieu use comconst_mod, only: kappa, cpp use comvert_mod, only: preff, pressure_exner, bp, ap USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h USE comgeom2_mod_h IMPLICIT NONE 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 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 CHARACTER(LEN=20),PARAMETER :: modname="guide_interp" write(*,*)trim(modname)//': interpolate nudging variables' ! ----------------------------------------------------------------- ! Calcul des niveaux de pression champs guidage ! ----------------------------------------------------------------- if (guide_modele) then do i=1,iip1 do j=1,jjp1 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=1,jjp1 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. write(*,*)trim(modname)//' : check vertical level order' write(*,*)trim(modname)//' LMDZ :' do l=1,llm write(*,*)trim(modname)//' PL(',l,')=',(ap(l)+ap(l+1))/2. & +psi(1,jjp1)*(bp(l)+bp(l+1))/2. enddo write(*,*)trim(modname)//' nudging file :' do l=1,nlevnc write(*,*)trim(modname)//' PL(',l,')=',plnc2(1,1,l) enddo write(*,*)trim(modname)//' invert ordering: invert_p=',invert_p if (guide_u) then do l=1,nlevnc write(*,*)trim(modname)//' U(',l,')=',unat2(1,1,l) enddo endif if (guide_T) then do l=1,nlevnc write(*,*)trim(modname)//' T(',l,')=',tnat2(1,1,l) enddo endif endif ! ----------------------------------------------------------------- ! Calcul niveaux pression modele ! ----------------------------------------------------------------- CALL pression( ip1jmp1, ap, bp, psi, p ) if (pressure_exner) then CALL exner_hyb(ip1jmp1,psi,p,pks,pk) else CALL exner_milieu(ip1jmp1,psi,p,pks,pk) endif ! .... Calcul de pls , pression au milieu des couches ,en Pascals unskap=1./kappa prefkap = preff ** kappa 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 ! calcul des pressions pour les grilles u et v 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,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=1,jjm 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=1,jjp1 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,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p) CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p) do l=1,llm do j=1,jjp1 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,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p) CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p) do l=1,llm do j=1,jjp1 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,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p) CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p) do l=1,llm do j=1,jjm 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,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p) CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p) do l=1,llm do j=1,jjp1 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*jjp1*llm,teta*pk/cpp,plsnc,qsat) qgui1=qgui1*qsat*0.01 !hum. rel. en % qgui2=qgui2*qsat*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) use comconst_mod, only: pi use serre_mod, only: clon, clat, grossismx, grossismy USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h USE comgeom2_mod_h implicit none ! 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 character(len=20),parameter :: modname="tau2alpha" CHARACTER (len = 80) :: abort_message 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 write(*,*)trim(modname)//' ATTENTION modele peu zoome' write(*,*)trim(modname)//' ATTENTION on prend une constante de guidage cste' gamma=0. else gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min) write(*,*)trim(modname)//' gamma=',gamma if (gamma.lt.1.e-5) then write(*,*)trim(modname)//' gamma =',gamma,'<1e-5' abort_message='stopped' CALL abort_gcm(modname,abort_message,1) endif gamma=log(0.5)/log(gamma) if (gamma4) then gamma=min(gamma,4.) endif write(*,*)trim(modname)//' 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 if (.not. guide_add) alpha = 1. - exp(- alpha) END SUBROUTINE tau2alpha !======================================================================= SUBROUTINE guide_read(timestep) use netcdf, only: NF90_GET_VAR, nf90_noerr USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h IMPLICIT NONE INTEGER, INTENT(IN) :: timestep LOGICAL, SAVE :: first=.TRUE. ! Identification fichiers et variables NetCDF: INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp INTEGER, SAVE :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps INTEGER :: ncidpl,varidpl,varidap,varidbp,dimid,lendim ! Variables auxiliaires NetCDF: INTEGER, DIMENSION(4) :: start,count INTEGER :: status,rcode CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: modname = 'guide_read' CHARACTER (len = 20) :: namedim ! ----------------------------------------------------------------- ! Premier appel: initialisation de la lecture des fichiers ! ----------------------------------------------------------------- if (first) then ncidpl=-99 write(*,*) trim(modname)//': opening nudging files ' ! Niveaux de pression si non constants if (guide_plevs.EQ.1) then write(*,*) trim(modname)//' Reading nudging on model levels' rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file apbp.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidpl, 'AP', varidap) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no AP variable in file apbp.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidpl, 'BP', varidbp) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no BP variable in file apbp.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap endif ! Pression si guidage sur niveaux P variables if (guide_plevs.EQ.2) then rcode = nf90_open('P.nc', nf90_nowrite, ncidp) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file P.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidp, 'PRES', varidp) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no PRES variable in file P.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp if (ncidpl.eq.-99) ncidpl=ncidp endif ! Vent zonal if (guide_u) then rcode = nf90_open('u.nc', nf90_nowrite, ncidu) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file u.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidu, 'UWND', varidu) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no UWND variable in file u.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu if (ncidpl.eq.-99) ncidpl=ncidu status=NF90_INQ_DIMID(ncidu, "LONU", dimid) status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim) IF (lendim .NE. iip1) THEN abort_message='dimension LONU different from iip1 in u.nc' CALL abort_gcm(modname,abort_message,1) ENDIF status=NF90_INQ_DIMID(ncidu, "LATU", dimid) status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim) IF (lendim .NE. jjp1) THEN abort_message='dimension LATU different from jjp1 in u.nc' CALL abort_gcm(modname,abort_message,1) ENDIF endif ! Vent meridien if (guide_v) then rcode = nf90_open('v.nc', nf90_nowrite, ncidv) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file v.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidv, 'VWND', varidv) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no VWND variable in file v.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv if (ncidpl.eq.-99) ncidpl=ncidv status=NF90_INQ_DIMID(ncidv, "LONV", dimid) status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim) IF (lendim .NE. iip1) THEN abort_message='dimension LONV different from iip1 in v.nc' CALL abort_gcm(modname,abort_message,1) ENDIF status=NF90_INQ_DIMID(ncidv, "LATV", dimid) status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim) IF (lendim .NE. jjm) THEN abort_message='dimension LATV different from jjm in v.nc' CALL abort_gcm(modname,abort_message,1) ENDIF endif ! Temperature if (guide_T) then rcode = nf90_open('T.nc', nf90_nowrite, ncidt) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file T.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidt, 'AIR', varidt) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no AIR variable in file T.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt if (ncidpl.eq.-99) ncidpl=ncidt status=NF90_INQ_DIMID(ncidt, "LONV", dimid) status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim) IF (lendim .NE. iip1) THEN abort_message='dimension LONV different from iip1 in T.nc' CALL abort_gcm(modname,abort_message,1) ENDIF status=NF90_INQ_DIMID(ncidt, "LATU", dimid) status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim) IF (lendim .NE. jjp1) THEN abort_message='dimension LATU different from jjp1 in T.nc' CALL abort_gcm(modname,abort_message,1) ENDIF endif ! Humidite if (guide_Q) then rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file hur.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidQ, 'RH', varidQ) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no RH variable in file hur.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ if (ncidpl.eq.-99) ncidpl=ncidQ status=NF90_INQ_DIMID(ncidQ, "LONV", dimid) status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim) IF (lendim .NE. iip1) THEN abort_message='dimension LONV different from iip1 in hur.nc' CALL abort_gcm(modname,abort_message,1) ENDIF status=NF90_INQ_DIMID(ncidQ, "LATU", dimid) status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim) IF (lendim .NE. jjp1) THEN abort_message='dimension LATU different from jjp1 in hur.nc' CALL abort_gcm(modname,abort_message,1) ENDIF endif ! Pression de surface if ((guide_P).OR.(guide_modele)) then rcode = nf90_open('ps.nc', nf90_nowrite, ncidps) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file ps.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidps, 'SP', varidps) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no SP variable in file ps.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps endif ! Coordonnee verticale if (guide_plevs.EQ.0) then rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl) IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl) write(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl endif ! Coefs ap, bp pour calcul de la pression aux differents niveaux if (guide_plevs.EQ.1) then status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc]) status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc]) ELSEIF (guide_plevs.EQ.0) THEN status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc]) !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous IF(convert_Pa) 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 ! Pression if (guide_plevs.EQ.2) then status=NF90_GET_VAR(ncidp,varidp,pnat2,start,count) IF (invert_y) THEN ! PRINT*,"Invertion impossible actuellement" ! CALL abort_gcm(modname,abort_message,1) CALL invert_lat(iip1,jjp1,nlevnc,pnat2) ENDIF endif ! Vent zonal if (guide_u) then status=NF90_GET_VAR(ncidu,varidu,unat2,start,count) IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,unat2) ENDIF endif ! Temperature if (guide_T) then status=NF90_GET_VAR(ncidt,varidt,tnat2,start,count) IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,tnat2) ENDIF endif ! Humidite if (guide_Q) then status=NF90_GET_VAR(ncidQ,varidQ,qnat2,start,count) IF (invert_y) THEN CALL invert_lat(iip1,jjp1,nlevnc,qnat2) ENDIF endif ! Vent meridien if (guide_v) then count(2)=jjm status=NF90_GET_VAR(ncidv,varidv,vnat2,start,count) 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 status=NF90_GET_VAR(ncidps,varidps,psnat2,start,count) IF (invert_y) THEN CALL invert_lat(iip1,jjp1,1,psnat2) ENDIF endif END SUBROUTINE guide_read !======================================================================= SUBROUTINE guide_read2D(timestep) use netcdf, only: nf90_get_var, nf90_noerr USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h IMPLICIT NONE INTEGER, INTENT(IN) :: timestep LOGICAL, SAVE :: first=.TRUE. ! Identification fichiers et variables NetCDF: INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp INTEGER, SAVE :: ncidQ,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 CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: modname = 'guide_read2D' ! ----------------------------------------------------------------- ! Premier appel: initialisation de la lecture des fichiers ! ----------------------------------------------------------------- if (first) then ncidpl=-99 write(*,*)trim(modname)//' : opening nudging files ' ! Ap et Bp si niveaux de pression hybrides if (guide_plevs.EQ.1) then write(*,*)trim(modname)//' Reading nudging on model levels' rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file apbp.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidpl, 'AP', varidap) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no AP variable in file apbp.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidpl, 'BP', varidbp) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no BP variable in file apbp.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*)trim(modname)//'ncidpl,varidap',ncidpl,varidap endif ! Pression if (guide_plevs.EQ.2) then rcode = nf90_open('P.nc', nf90_nowrite, ncidp) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file P.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidp, 'PRES', varidp) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no PRES variable in file P.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*)trim(modname)//' ncidp,varidp',ncidp,varidp if (ncidpl.eq.-99) ncidpl=ncidp endif ! Vent zonal if (guide_u) then rcode = nf90_open('u.nc', nf90_nowrite, ncidu) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file u.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidu, 'UWND', varidu) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no UWND variable in file u.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*)trim(modname)//' 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) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file v.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidv, 'VWND', varidv) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no VWND variable in file v.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*)trim(modname)//' ncidv,varidv',ncidv,varidv if (ncidpl.eq.-99) ncidpl=ncidv endif ! Temperature if (guide_T) then rcode = nf90_open('T.nc', nf90_nowrite, ncidt) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file T.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidt, 'AIR', varidt) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no AIR variable in file T.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*)trim(modname)//' ncidT,varidT',ncidt,varidt if (ncidpl.eq.-99) ncidpl=ncidt endif ! Humidite if (guide_Q) then rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file hur.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidQ, 'RH', varidQ) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no RH,variable in file hur.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*)trim(modname)//' 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) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no file ps.nc' CALL abort_gcm(modname,abort_message,1) ENDIF rcode = nf90_inq_varid(ncidps, 'SP', varidps) IF (rcode.NE.NF90_NOERR) THEN abort_message='Nudging: error -> no SP variable in file ps.nc' CALL abort_gcm(modname,abort_message,1) ENDIF write(*,*)trim(modname)//' ncidps,varidps',ncidps,varidps endif ! Coordonnee verticale if (guide_plevs.EQ.0) then rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl) IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl) write(*,*)trim(modname)//' ncidpl,varidpl',ncidpl,varidpl endif ! Coefs ap, bp pour calcul de la pression aux differents niveaux if (guide_plevs.EQ.1) then status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc]) status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc]) elseif (guide_plevs.EQ.0) THEN status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc]) 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 ! Pression if (guide_plevs.EQ.2) then status=NF90_GET_VAR(ncidp,varidp,zu,start,count) DO i=1,iip1 pnat2(i,:,:)=zu(:,:) ENDDO IF (invert_y) THEN ! PRINT*,"Invertion impossible actuellement" ! CALL abort_gcm(modname,abort_message,1) CALL invert_lat(iip1,jjp1,nlevnc,pnat2) ENDIF endif ! Vent zonal if (guide_u) then status=NF90_GET_VAR(ncidu,varidu,zu,start,count) 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 status=NF90_GET_VAR(ncidt,varidt,zu,start,count) 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 status=NF90_GET_VAR(ncidQ,varidQ,zu,start,count) 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 status=NF90_GET_VAR(ncidv,varidv,zv,start,count) 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_plevs.EQ.1)) then start(3)=timestep start(4)=0 count(2)=jjp1 count(3)=1 count(4)=0 status=NF90_GET_VAR(ncidps,varidps,zu(:,1),start,count) 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 comconst_mod, ONLY: pi USE comvert_mod, ONLY: presnivs use netcdf95, only: nf95_def_var, nf95_put_var use netcdf, only: nf90_float, nf90_def_var, nf90_put_var USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h USE comgeom2_mod_h IMPLICIT NONE ! 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 :: vid_au,vid_av, varid_alpha_t, varid_alpha_q INTEGER, DIMENSION (3) :: dim3 INTEGER, DIMENSION (4) :: dim4,count,start INTEGER :: ierr, varid,l REAL, DIMENSION (iip1,hsize,vsize) :: field2 CHARACTER(LEN=20),PARAMETER :: modname="guide_out" write(*,*)trim(modname)//': output timestep',timestep,'var ',varname IF (timestep.EQ.0) THEN ! ---------------------------------------------- ! initialisation fichier de sortie ! ---------------------------------------------- ! Ouverture du fichier ierr=nf90_create("guide_ins.nc",IOR(nf90_clobber,nf90_64bit_offset),nid) ! Definition des dimensions ierr=nf90_def_dim(nid,"LONU",iip1,id_lonu) ierr=nf90_def_dim(nid,"LONV",iip1,id_lonv) ierr=nf90_def_dim(nid,"LATU",jjp1,id_latu) ierr=nf90_def_dim(nid,"LATV",jjm,id_latv) ierr=nf90_def_dim(nid,"LEVEL",llm,id_lev) ierr=nf90_def_dim(nid,"TIME",nf90_unlimited,id_tim) ! Creation des variables dimensions ierr=NF90_DEF_VAR(nid,"LONU",NF90_FLOAT,id_lonu,vid_lonu) ierr=NF90_DEF_VAR(nid,"LONV",NF90_FLOAT,id_lonv,vid_lonv) ierr=NF90_DEF_VAR(nid,"LATU",NF90_FLOAT,id_latu,vid_latu) ierr=NF90_DEF_VAR(nid,"LATV",NF90_FLOAT,id_latv,vid_latv) ierr=NF90_DEF_VAR(nid,"LEVEL",NF90_FLOAT,id_lev,vid_lev) ierr=NF90_DEF_VAR(nid,"cu",NF90_FLOAT,(/id_lonu,id_latu/),vid_cu) ierr=NF90_DEF_VAR(nid,"cv",NF90_FLOAT,(/id_lonv,id_latv/),vid_cv) ierr=NF90_DEF_VAR(nid,"au",NF90_FLOAT,(/id_lonu,id_latu/),vid_au) ierr=NF90_DEF_VAR(nid,"av",NF90_FLOAT,(/id_lonv,id_latv/),vid_av) call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), & varid_alpha_t) call nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), & varid_alpha_q) ierr=nf90_enddef(nid) ! Enregistrement des variables dimensions ierr = nf90_put_var(nid, vid_lonu, rlonu * 180. / pi) ierr = nf90_put_var(nid, vid_lonv, rlonv * 180. / pi) ierr = nf90_put_var(nid, vid_latu, rlatu * 180. / pi) ierr = nf90_put_var(nid, vid_latv, rlatv * 180. / pi) ierr = nf90_put_var(nid, vid_lev, presnivs) ierr = nf90_put_var(nid, vid_cu, cu) ierr = nf90_put_var(nid, vid_cv, cv) ierr = nf90_put_var(nid, vid_au, alpha_u) ierr = nf90_put_var(nid, vid_av, alpha_v) call nf95_put_var(nid, varid_alpha_t, alpha_t) call nf95_put_var(nid, varid_alpha_q, alpha_q) ! -------------------------------------------------------------------- ! Cr�ation des variables sauvegard�es ! -------------------------------------------------------------------- ierr = nf90_redef(nid) ! Pressure (GCM) dim4=(/id_lonv,id_latu,id_lev,id_tim/) ierr = NF90_DEF_VAR(nid,"SP",NF90_FLOAT,dim4,varid) ! Surface pressure (guidage) IF (guide_P) THEN dim3=(/id_lonv,id_latu,id_tim/) ierr = NF90_DEF_VAR(nid,"ps",NF90_FLOAT,dim3,varid) ENDIF ! Zonal wind IF (guide_u) THEN dim4=(/id_lonu,id_latu,id_lev,id_tim/) ierr = NF90_DEF_VAR(nid,"u",NF90_FLOAT,dim4,varid) ierr = NF90_DEF_VAR(nid,"ua",NF90_FLOAT,dim4,varid) ierr = NF90_DEF_VAR(nid,"ucov",NF90_FLOAT,dim4,varid) ENDIF ! Merid. wind IF (guide_v) THEN dim4=(/id_lonv,id_latv,id_lev,id_tim/) ierr = NF90_DEF_VAR(nid,"v",NF90_FLOAT,dim4,varid) ierr = NF90_DEF_VAR(nid,"va",NF90_FLOAT,dim4,varid) ierr = NF90_DEF_VAR(nid,"vcov",NF90_FLOAT,dim4,varid) ENDIF ! Pot. Temperature IF (guide_T) THEN dim4=(/id_lonv,id_latu,id_lev,id_tim/) ierr = NF90_DEF_VAR(nid,"teta",NF90_FLOAT,dim4,varid) ENDIF ! Specific Humidity IF (guide_Q) THEN dim4=(/id_lonv,id_latu,id_lev,id_tim/) ierr = NF90_DEF_VAR(nid,"q",NF90_FLOAT,dim4,varid) ENDIF ierr = nf90_enddef(nid) ierr = nf90_close(nid) ENDIF ! timestep=0 ! -------------------------------------------------------------------- ! Enregistrement du champ ! -------------------------------------------------------------------- ierr=nf90_open("guide_ins.nc",nf90_write,nid) IF (varname=="SP") timestep=timestep+1 ierr = nf90_inq_varid(nid,varname,varid) SELECT CASE (varname) CASE ("SP","ps") start=(/1,1,1,timestep/) count=(/iip1,jjp1,llm,1/) CASE ("v","va","vcov") start=(/1,1,1,timestep/) count=(/iip1,jjm,llm,1/) CASE DEFAULT start=(/1,1,1,timestep/) count=(/iip1,jjp1,llm,1/) END SELECT SELECT CASE (varname) CASE("u","ua") DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO field2(:,1,:)=0. ; field2(:,jjp1,:)=0. CASE("v","va") DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO CASE DEFAULT field2=field END SELECT ierr = nf90_put_var(nid, varid, field2, start, count) ierr = nf90_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_mod