MODULE guide_loc_mod !======================================================================= ! Auteur: F.Hourdin ! F. Codron 01/09 !======================================================================= USE getparam, ONLY: ini_getparam, fin_getparam, getpar USE Write_Field_loc use netcdf, ONLY: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, & nf90_inq_dimid, nf90_inquire_dimension, nf90_inq_dimid, & nf90_inquire_dimension, nf90_enddef, nf90_def_dim, nf90_put_var, nf90_noerr, nf90_close, nf90_inq_varid, & nf90_redef, nf90_write, nf90_unlimited, nf90_float, nf90_clobber, nf90_64bit_offset, nf90_float, & nf90_create, nf90_def_var, nf90_open USE parallel_lmdz USE pres2lev_mod, ONLY: pres2lev 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 INTEGER,SAVE,PRIVATE :: ijbu,ijbv,ijeu,ijev !,ijnu,ijnv INTEGER,SAVE,PRIVATE :: jjbu,jjbv,jjeu,jjev,jjnu,jjnv CONTAINS !======================================================================= SUBROUTINE guide_init USE control_mod, ONLY: day_step USE serre_mod, ONLY: grossismx IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" 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>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==0) THEN iguide_int=1 ELSEIF (iguide_int>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==1) THEN if (ncidpl==-99) THEN rcod=nf90_open('apbp.nc',nf90_nowrite, ncidpl) if (rcod/=nf90_noerr) THEN abort_message=' Nudging error -> no file apbp.nc' CALL abort_gcm(modname,abort_message,1) endif endif elseif (guide_plevs==2) THEN if (ncidpl==-99) THEN rcod=nf90_open('P.nc',nf90_nowrite,ncidpl) if (rcod/=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==-99) THEN rcod=nf90_open('u.nc',nf90_nowrite,ncidpl) if (rcod/=nf90_noerr) THEN abort_message=' Nudging error -> no file u.nc' CALL abort_gcm(modname,abort_message,1) endif endif elseif (guide_v) THEN if (ncidpl==-99) THEN rcod=nf90_open('v.nc',nf90_nowrite,ncidpl) if (rcod/=nf90_noerr) THEN abort_message=' Nudging error -> no file v.nc' CALL abort_gcm(modname,abort_message,1) endif endif elseif (guide_T) THEN if (ncidpl==-99) THEN rcod=nf90_open('T.nc',nf90_nowrite,ncidpl) if (rcod/=nf90_noerr) THEN abort_message=' Nudging error -> no file T.nc' CALL abort_gcm(modname,abort_message,1) endif endif elseif (guide_Q) THEN if (ncidpl==-99) THEN rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl) if (rcod/=nf90_noerr) THEN abort_message=' Nudging error -> no file hur.nc' CALL abort_gcm(modname,abort_message,1) endif endif endif error=nf90_inq_dimid(ncidpl,'LEVEL',rid) IF (error/=nf90_noerr) error=nf90_inq_dimid(ncidpl,'PRESSURE',rid) IF (error/=nf90_noerr) THEN abort_message='Nudging: error reading pressure levels' CALL abort_gcm(modname,abort_message,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(ijb_u:ije_u), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_v(ijb_v:ije_v), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_T(ijb_u:ije_u), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_Q(ijb_u:ije_u), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(alpha_P(ijb_u:ije_u), 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,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(ugui1(ijb_u:ije_u,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(unat2(iip1,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(ugui2(ijb_u:ije_u,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,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(tgui1(ijb_u:ije_u,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(tnat2(iip1,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(tgui2(ijb_u:ije_u,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,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(qgui1(ijb_u:ije_u,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(qnat2(iip1,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(qgui2(ijb_u:ije_u,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,jjb_v:jje_v,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(vgui1(ijb_v:ije_v,llm), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(vnat2(iip1,jjb_v:jje_v,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(vgui2(ijb_v:ije_v,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==2) THEN ALLOCATE(pnat1(iip1,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(pnat2(iip1,jjb_u:jje_u,nlevnc), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) pnat1=0.;pnat2=0.; ENDIF IF (guide_P.OR.guide_plevs==1) THEN ALLOCATE(psnat1(iip1,jjb_u:jje_u), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(psnat2(iip1,jjb_u:jje_u), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) psnat1=0.;psnat2=0.; ENDIF IF (guide_P) THEN ALLOCATE(psgui2(ijb_u:ije_u), stat = error) IF (error /= 0) CALL abort_gcm(modname,abort_message,1) ALLOCATE(psgui1(ijb_u:ije_u), 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==2) pnat1=pnat2 IF (guide_P.OR.guide_plevs==1) psnat1=psnat2 END SUBROUTINE guide_init !======================================================================= SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps) use exner_hyb_loc_m, ONLY: exner_hyb_loc use exner_milieu_loc_m, ONLY: exner_milieu_loc USE parallel_lmdz USE control_mod USE write_field_loc USE comconst_mod, ONLY: cpp, daysec, dtvr, kappa USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" ! Variables entree INTEGER, INTENT(IN) :: itau !pas de temps REAL, DIMENSION (ijb_u:ije_u,llm), INTENT(INOUT) :: ucov,teta,q,masse REAL, DIMENSION (ijb_v:ije_v,llm), INTENT(INOUT) :: vcov REAL, DIMENSION (ijb_u:ije_u), INTENT(INOUT) :: ps ! Variables locales LOGICAL, SAVE :: first=.TRUE. !$OMP THREADPRIVATE(first) LOGICAL :: f_out ! sortie guidage REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addu ! var aux: champ de guidage REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addv ! var aux: champ de guidage ! Variables pour fonction Exner (P milieu couche) REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: pk REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: pks REAL :: unskap REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: p ! besoin si guide_P ! Compteurs temps: INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage !$OMP THREADPRIVATE(step_rea,count_no_rea,itau_test) REAL :: ditau, dday_step REAL :: tau,reste ! position entre 2 etats de guidage REAL, SAVE :: factt ! pas de temps en fraction de jour !$OMP THREADPRIVATE(factt) INTEGER :: i,j,l CHARACTER(LEN=20) :: modname="guide_main" !$OMP MASTER ijbu=ij_begin ; ijeu=ij_end jjbu=jj_begin ; jjeu=jj_end ; jjnu=jjeu-jjbu+1 ijbv=ij_begin ; ijev=ij_end jjbv=jj_begin ; jjev=jj_end ; jjnv=jjev-jjbv+1 IF (pole_sud) THEN ijeu=ij_end-iip1 ijev=ij_end-iip1 jjev=jj_end-1 jjnv=jjev-jjbv+1 ENDIF IF (pole_nord) THEN ijbu=ij_begin+iip1 ijbv=ij_begin ENDIF !$OMP END MASTER !$OMP BARRIER ! 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. !$OMP MASTER ALLOCATE(f_addu(ijb_u:ije_u,llm) ) ALLOCATE(f_addv(ijb_v:ije_v,llm) ) ALLOCATE(pk(iip1,jjb_u:jje_u,llm) ) ALLOCATE(pks(iip1,jjb_u:jje_u) ) ALLOCATE(p(ijb_u:ije_u,llmp1) ) CALL guide_init !$OMP END MASTER !$OMP BARRIER itau_test=1001 step_rea=1 count_no_rea=0 ! Calcul des constantes de rappel factt=dtvr*iperiod/daysec !$OMP MASTER CALL tau2alpha(3, iip1, jjb_v, jje_v, factt, tau_min_v, tau_max_v, alpha_v) CALL tau2alpha(2, iip1, jjb_u, jje_u, factt, tau_min_u, tau_max_u, alpha_u) CALL tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_T, tau_max_T, alpha_T) CALL tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_P, tau_max_P, alpha_P) CALL tau2alpha(1, iip1, jjb_u, jje_u, 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 !$OMP END MASTER !$OMP BARRIER ! ini_anal: etat initial egal au guidage IF (ini_anal) THEN CALL guide_interp(ps,teta) !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm IF (guide_u) ucov(ijbu:ijeu,l)=ugui2(ijbu:ijeu,l) IF (guide_v) vcov(ijbv:ijev,l)=ugui2(ijbv:ijev,l) IF (guide_T) teta(ijbu:ijeu,l)=tgui2(ijbu:ijeu,l) IF (guide_Q) q(ijbu:ijeu,l)=qgui2(ijbu:ijeu,l) ENDDO IF (guide_P) THEN !$OMP MASTER ps(ijbu:ijeu)=psgui2(ijbu:ijeu) !$OMP END MASTER !$OMP BARRIER CALL pression_loc(ijnb_u,ap,bp,ps,p) CALL massdair_loc(p,masse) !$OMP BARRIER ENDIF RETURN ENDIF ENDIF !first !----------------------------------------------------------------------- ! Lecture des fichiers de guidage ? !----------------------------------------------------------------------- IF (iguide_read/=0) THEN ditau=real(itau) dday_step=real(day_step) IF (iguide_read<0) THEN tau=ditau/dday_step/REAL(iguide_read) ELSE tau=REAL(iguide_read)*ditau/dday_step ENDIF reste=tau-AINT(tau) IF (reste==0.) THEN IF (itau_test==itau) THEN WRITE(*,*)trim(modname)//' second pass in advreel at itau=',& itau CALL abort_gcm("guide_loc_lod","stopped",1) ELSE !$OMP MASTER IF (guide_v) vnat1(:,jjbv:jjev,:)=vnat2(:,jjbv:jjev,:) IF (guide_u) unat1(:,jjbu:jjeu,:)=unat2(:,jjbu:jjeu,:) IF (guide_T) tnat1(:,jjbu:jjeu,:)=tnat2(:,jjbu:jjeu,:) IF (guide_Q) qnat1(:,jjbu:jjeu,:)=qnat2(:,jjbu:jjeu,:) IF (guide_plevs==2) pnat1(:,jjbu:jjeu,:)=pnat2(:,jjbu:jjeu,:) IF (guide_P.OR.guide_plevs==1) psnat1(:,jjbu:jjeu)=psnat2(:,jjbu:jjeu) !$OMP END MASTER !$OMP BARRIER step_rea=step_rea+1 itau_test=itau if (is_master) THEN WRITE(*,*)trim(modname)//' Reading nudging files, step ',& step_rea,'after ',count_no_rea,' skips' endif IF (guide_2D) THEN !$OMP MASTER CALL guide_read2D(step_rea) !$OMP END MASTER !$OMP BARRIER ELSE !$OMP MASTER CALL guide_read(step_rea) !$OMP END MASTER !$OMP BARRIER 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)==0) THEN CALL guide_interp(ps,teta) ENDIF ! Repartition entre 2 etats de guidage IF (iguide_read/=0) THEN tau=reste ELSE tau=1. ENDIF ! CALL WriteField_u('ucov_guide',ucov) ! CALL WriteField_v('vcov_guide',vcov) ! CALL WriteField_u('teta_guide',teta) ! CALL WriteField_u('masse_guide',masse) !----------------------------------------------------------------------- ! Ajout des champs de guidage !----------------------------------------------------------------------- ! Sauvegarde du guidage? f_out=((MOD(itau,iguide_sav)==0).AND.guide_sav) IF (f_out) THEN !$OMP BARRIER CALL pression_loc(ijnb_u,ap,bp,ps,p) !$OMP BARRIER if (pressure_exner) THEN CALL exner_hyb_loc( ijnb_u, ps, p, pks, pk) else CALL exner_milieu_loc( ijnb_u, ps, p, pks, pk ) endif !$OMP BARRIER unskap=1./kappa !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO j=jjbu,jjeu DO i =1, iip1 p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap ENDDO ENDDO ENDDO CALL guide_out("SP",jjp1,llm,p(ijb_u:ije_u,1:llm),1.) ENDIF if (guide_u) THEN if (guide_add) THEN !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm f_addu(ijbu:ijeu,l)=(1.-tau)*ugui1(ijbu:ijeu,l)+tau*ugui2(ijbu:ijeu,l) ENDDO else !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm f_addu(ijbu:ijeu,l)=(1.-tau)*ugui1(ijbu:ijeu,l)+tau*ugui2(ijbu:ijeu,l)-ucov(ijbu:ijeu,l) ENDDO endif ! CALL WriteField_u('f_addu',f_addu) if (guide_zon) CALL guide_zonave_u(1,llm,f_addu) CALL guide_addfield_u(llm,f_addu,alpha_u) IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1(ijb_u:ije_u,:)+tau*ugui2(ijb_u:ije_u,:),factt) IF (f_out) CALL guide_out("u",jjp1,llm,ucov(ijb_u:ije_u,:),factt) IF (f_out) THEN ! Ehouarn: fill the gaps adequately... IF (ijbu>ijb_u) f_addu(ijb_u:ijbu-1,:)=0 IF (ijeuijb_v) f_addv(ijb_v:ijbv-1,:)=0 IF (ijev a verifier ! ym DO j=jjbv,jjev DO j=jjbu,jjeu 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=jjbu,jjeu DO i=1,iip1 ij=(j-1)*iip1+i field(ij,l)=fieldm(j,l) ENDDO ENDDO ENDDO END SUBROUTINE guide_zonave_u SUBROUTINE guide_zonave_v(typ,hsize,vsize,field) USE comconst_mod, ONLY: pi IMPLICIT NONE INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "comgeom.h" ! input/output variables INTEGER, INTENT(IN) :: typ INTEGER, INTENT(IN) :: vsize INTEGER, INTENT(IN) :: hsize REAL, DIMENSION(ijb_v:ije_v,vsize), INTENT(INOUT) :: field ! Local variables LOGICAL, SAVE :: first=.TRUE. !$OMP THREADPRIVATE(first) INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain !$OMP THREADPRIVATE(imin, imax) INTEGER :: i,j,l,ij REAL, DIMENSION (iip1) :: lond ! longitude in Deg. REAL, DIMENSION (jjb_v:jjev,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) 0) THEN !$OMP MASTER DEALLOCATE(field_glo) !$OMP END MASTER !$OMP BARRIER RETURN ENDIF !$OMP MASTER IF (timestep==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,zu) ierr = nf90_put_var(nid,vid_av,zv) CALL nf95_put_var(nid, varid_alpha_t, zt) CALL nf95_put_var(nid, varid_alpha_q, zq) ! -------------------------------------------------------------------- ! 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 !$OMP END MASTER !$OMP BARRIER SELECT CASE (varname) CASE("u","ua") !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm field_glo(:,2:jjm,l)=field_glo(:,2:jjm,l)/cu(:,2:jjm) field_glo(:,1,l)=0. ; field_glo(:,jjp1,l)=0. ENDDO CASE("v","va") !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm field_glo(:,:,l)=field_glo(:,:,l)/cv(:,:) ENDDO END SELECT ! if (varname=="ua") THEN ! CALL dump2d(iip1,jjp1,field_glo,'ua gui1 1ere couche ') ! CALL dump2d(iip1,jjp1,field_glo(:,:,llm),'ua gui1 llm ') ! endif !$OMP MASTER ierr = nf90_put_var(nid,varid,field_glo,start,count) ierr = nf90_close(nid) DEALLOCATE(field_glo) !$OMP END MASTER !$OMP BARRIER 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))>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 END SUBROUTINE correctbid !==================================================================== ! Ascii debug output. Could be reactivated !==================================================================== SUBROUTINE dump2du(var,varname) use parallel_lmdz use mod_hallo IMPLICIT NONE include 'dimensions.h' include 'paramet.h' CHARACTER (len=*) :: varname real, dimension(ijb_u:ije_u) :: var real, dimension(ip1jmp1) :: var_glob RETURN CALL barrier CALL Gather_field_u(var,var_glob,1) CALL barrier if (mpi_rank==0) THEN CALL dump2d(iip1,jjp1,var_glob,varname) endif CALL barrier END SUBROUTINE dump2du !==================================================================== ! Ascii debug output. Could be reactivated !==================================================================== SUBROUTINE dumpall IMPLICIT NONE include "dimensions.h" include "paramet.h" include "comgeom.h" CALL barrier CALL dump2du(alpha_u(ijb_u:ije_u),' alpha_u couche 1') CALL dump2du(unat2(:,jjbu:jjeu,nlevnc),' unat2 couche nlevnc') CALL dump2du(ugui1(ijb_u:ije_u,1)*sqrt(unscu2(ijb_u:ije_u)),' ugui1 couche 1') END SUBROUTINE dumpall !=========================================================================== END MODULE guide_loc_mod