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 lmdz_pres2lev, 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 USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet 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>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 USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! 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 USE lmdz_comgeom USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! 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 USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE 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 USE lmdz_comgeom USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE 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