MODULE lmdz_old_1dconv PRIVATE ! -- We'd love to put IMPLICIT NONE; here... PUBLIC get_uvd, copie, get_uvd2, rdgrads, spaces CONTAINS subroutine get_uvd(itap, dtime, file_forctl, file_fordat, & & ht, hq, hw, hu, hv, hthturb, hqturb, & & Ts, imp_fcg, ts_fcg, Tp_fcg, Turb_fcg) USE lmdz_yomcst IMPLICIT NONE !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de ! pouvoir calculer la convergence et le cisaillement dans la physiq !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER klev REAL play(100) !pression en Pa au milieu de chaque couche GCM INTEGER JM(100) !pression en Pa au milieu de chaque couche GCM REAL coef1(100) !coefficient d interpolation REAL coef2(100) !coefficient d interpolation INTEGER nblvlm !nombre de niveau de pression du mesoNH REAL playm(100) !pression en Pa au milieu de chaque couche Meso-NH REAL hplaym(100) !pression en hPa milieux des couches Meso-NH INTEGER i, j, k, ll, in CHARACTER*80 file_forctl, file_fordat COMMON/com1_phys_gcss/play, coef1, coef2, JM, klev COMMON/com2_phys_gcss/playm, hplaym, nblvlm !====================================================================== ! methode: on va chercher les donnees du mesoNH de meteo france, on y ! a acces a tout pas detemps grace a la routine rdgrads qui ! est une boucle lisant dans ces fichiers. ! Puis on interpole ces donnes sur les 11 niveaux du gcm et ! et sur les pas de temps de ce meme gcm !---------------------------------------------------------------------- ! input: ! pasmax :nombre de pas de temps maximum du mesoNH ! dt :pas de temps du meso_NH (en secondes) !---------------------------------------------------------------------- INTEGER pasmax, dt save pasmax, dt !---------------------------------------------------------------------- ! arguments: ! itap :compteur de la physique(le nombre de ces pas est ! fixe dans la subroutine calcul_ini_gcm de interpo ! -lation ! dtime :pas detemps du gcm (en secondes) ! ht :convergence horizontale de temperature(K/s) ! hq : " " d humidite (kg/kg/s) ! hw :vitesse verticale moyenne (m/s**2) ! hu :convergence horizontale d impulsion le long de x ! (kg/(m^2 s^2) ! hv : idem le long de y. ! Ts : Temperature de surface (K) ! imp_fcg: var. LOGICAL .EQ. T si forcage en impulsion ! ts_fcg: var. LOGICAL .EQ. T si forcage en Ts present dans fichier ! Tp_fcg: var. LOGICAL .EQ. T si forcage donne en Temp potentielle ! Turb_fcg: var. LOGICAL .EQ. T si forcage turbulent present dans fichier !---------------------------------------------------------------------- INTEGER itap REAL dtime REAL ht(:) REAL hq(:) REAL hu(:) REAL hv(:) REAL hw(:) REAL hthturb(:) REAL hqturb(:) REAL Ts, Ts_subr LOGICAL imp_fcg LOGICAL ts_fcg LOGICAL Tp_fcg LOGICAL Turb_fcg !---------------------------------------------------------------------- ! Variables internes de get_uvd (note : l interpolation temporelle ! est faite entre les pas de temps before et after, sur les variables ! definies sur la grille du SCM; on atteint exactement les valeurs Meso ! aux milieux des pas de temps Meso) ! time0 :date initiale en secondes ! time :temps associe a chaque pas du SCM ! pas :numero du pas du meso_NH (on lit en pas : le premier pas ! des donnees est duplique) ! pasprev :numero du pas de lecture precedent ! htaft :advection horizontale de temp. au pas de temps after ! hqaft : " " d humidite " ! hwaft :vitesse verticalle moyenne au pas de temps after ! huaft,hvaft :advection horizontale d impulsion au pas de temps after ! tsaft : surface temperature 'after time step' ! htbef :idem htaft, mais pour le pas de temps before ! hqbef :voir hqaft ! hwbef :voir hwaft ! hubef,hvbef : idem huaft,hvaft, mais pour before ! tsbef : surface temperature 'before time step' !---------------------------------------------------------------------- INTEGER time0, pas, pasprev save time0, pas, pasprev REAL time REAL htaft(100), hqaft(100), hwaft(100), huaft(100), hvaft(100) REAL hthturbaft(100), hqturbaft(100) REAL Tsaft save htaft, hqaft, hwaft, huaft, hvaft, hthturbaft, hqturbaft REAL htbef(100), hqbef(100), hwbef(100), hubef(100), hvbef(100) REAL hthturbbef(100), hqturbbef(100) REAL Tsbef save htbef, hqbef, hwbef, hubef, hvbef, hthturbbef, hqturbbef REAL timeaft, timebef save timeaft, timebef INTEGER temps CHARACTER*4 string !---------------------------------------------------------------------- ! variables arguments de la subroutine rdgrads !--------------------------------------------------------------------- INTEGER icompt, icomp1 !compteurs de rdgrads REAL z(100) ! altitude (grille Meso) REAL ht_mes(100) !convergence horizontale de temperature !-(grille Meso) REAL hq_mes(100) !convergence horizontale d humidite !(grille Meso) REAL hw_mes(100) !vitesse verticale moyenne !(grille Meso) REAL hu_mes(100), hv_mes(100) !convergence horizontale d impulsion !(grille Meso) REAL hthturb_mes(100) !tendance horizontale de T_pot, due aux !flux turbulents REAL hqturb_mes(100) !tendance horizontale d humidite, due aux !flux turbulents !--------------------------------------------------------------------- ! variable argument de la subroutine copie !--------------------------------------------------------------------- ! SB real pplay(100) !pression en milieu de couche du gcm ! SB !argument de la physique !--------------------------------------------------------------------- ! variables destinees a la lecture du pas de temps du fichier de donnees !--------------------------------------------------------------------- CHARACTER*80 aaa, atemps, apasmax INTEGER nch, imn, ipa !--------------------------------------------------------------------- PRINT*, 'le pas itap est:', itap !*** on determine le pas du meso_NH correspondant au nouvel itap *** !*** pour aller chercher les champs dans rdgrads *** time = time0 + itap * dtime !c temps=int(time/dt+1) !c pas=min(temps,pasmax) temps = 1 + int((dt + 2 * time) / (2 * dt)) pas = min(temps, pasmax - 1) PRINT*, 'le pas Meso est:', pas !=================================================================== !*** on remplit les champs before avec les champs after du pas *** !*** precedent en format gcm *** IF(pas>pasprev)THEN do i = 1, klev htbef(i) = htaft(i) hqbef(i) = hqaft(i) hwbef(i) = hwaft(i) hubef(i) = huaft(i) hvbef(i) = hvaft(i) hThTurbbef(i) = hThTurbaft(i) hqTurbbef(i) = hqTurbaft(i) enddo tsbef = tsaft timebef = pasprev * dt timeaft = timebef + dt icomp1 = nblvlm * 4 IF (ts_fcg) icomp1 = icomp1 + 1 IF (imp_fcg) icomp1 = icomp1 + nblvlm * 2 IF (Turb_fcg) icomp1 = icomp1 + nblvlm * 2 icompt = icomp1 * pas print *, 'imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt' print *, imp_fcg, ts_fcg, Turb_fcg, pas, nblvlm, icompt PRINT*, 'le pas pas est:', pas !*** on va chercher les nouveaux champs after dans toga.dat *** !*** champs en format meso_NH *** open(99, FILE = file_fordat, FORM = 'UNFORMATTED', & & ACCESS = 'DIRECT', RECL = 8) CALL rdgrads(99, icompt, nblvlm, z, ht_mes, hq_mes, hw_mes & &, hu_mes, hv_mes, hthturb_mes, hqturb_mes & &, ts_fcg, ts_subr, imp_fcg, Turb_fcg) IF(Tp_fcg) THEN ! (le forcage est donne en temperature potentielle) do i = 1, nblvlm ht_mes(i) = ht_mes(i) * (hplaym(i) / 1000.)**rkappa enddo endif ! Tp_fcg IF(Turb_fcg) THEN do i = 1, nblvlm hThTurb_mes(i) = hThTurb_mes(i) * (hplaym(i) / 1000.)**rkappa enddo endif ! Turb_fcg PRINT*, 'ht_mes ', (ht_mes(i), i = 1, nblvlm) PRINT*, 'hq_mes ', (hq_mes(i), i = 1, nblvlm) PRINT*, 'hw_mes ', (hw_mes(i), i = 1, nblvlm) IF(imp_fcg) THEN PRINT*, 'hu_mes ', (hu_mes(i), i = 1, nblvlm) PRINT*, 'hv_mes ', (hv_mes(i), i = 1, nblvlm) endif IF(Turb_fcg) THEN PRINT*, 'hThTurb_mes ', (hThTurb_mes(i), i = 1, nblvlm) PRINT*, 'hqTurb_mes ', (hqTurb_mes(i), i = 1, nblvlm) endif IF (ts_fcg) PRINT*, 'ts_subr', ts_subr !*** on interpole les champs meso_NH sur les niveaux de pression*** !*** gcm . on obtient le nouveau champ after *** do k = 1, klev IF (JM(k) == 0) THEN htaft(k) = ht_mes(jm(k) + 1) hqaft(k) = hq_mes(jm(k) + 1) hwaft(k) = hw_mes(jm(k) + 1) IF(imp_fcg) THEN huaft(k) = hu_mes(jm(k) + 1) hvaft(k) = hv_mes(jm(k) + 1) endif ! imp_fcg IF(Turb_fcg) THEN hThTurbaft(k) = hThTurb_mes(jm(k) + 1) hqTurbaft(k) = hqTurb_mes(jm(k) + 1) endif ! Turb_fcg else ! JM(k) .EQ. 0 htaft(k) = coef1(k) * ht_mes(jm(k)) + coef2(k) * ht_mes(jm(k) + 1) hqaft(k) = coef1(k) * hq_mes(jm(k)) + coef2(k) * hq_mes(jm(k) + 1) hwaft(k) = coef1(k) * hw_mes(jm(k)) + coef2(k) * hw_mes(jm(k) + 1) IF(imp_fcg) THEN huaft(k) = coef1(k) * hu_mes(jm(k)) + coef2(k) * hu_mes(jm(k) + 1) hvaft(k) = coef1(k) * hv_mes(jm(k)) + coef2(k) * hv_mes(jm(k) + 1) endif ! imp_fcg IF(Turb_fcg) THEN hThTurbaft(k) = coef1(k) * hThTurb_mes(jm(k)) & & + coef2(k) * hThTurb_mes(jm(k) + 1) hqTurbaft(k) = coef1(k) * hqTurb_mes(jm(k)) & & + coef2(k) * hqTurb_mes(jm(k) + 1) endif ! Turb_fcg endif ! JM(k) .EQ. 0 enddo tsaft = ts_subr pasprev = pas else ! pas.gt.pasprev PRINT*, 'timebef est:', timebef endif ! pas.gt.pasprev fin du bloc relatif au passage au pas !de temps (meso) suivant !*** si on atteint le pas max des donnees experimentales ,on *** !*** on conserve les derniers champs calcules *** IF(temps>=pasmax)THEN do ll = 1, klev ht(ll) = htaft(ll) hq(ll) = hqaft(ll) hw(ll) = hwaft(ll) hu(ll) = huaft(ll) hv(ll) = hvaft(ll) hThTurb(ll) = hThTurbaft(ll) hqTurb(ll) = hqTurbaft(ll) enddo ts_subr = tsaft else ! temps.ge.pasmax !*** on interpole sur les pas de temps de 10mn du gcm a partir *** !** des pas de temps de 1h du meso_NH *** do j = 1, klev ht(j) = ((timeaft - time) * htbef(j) + (time - timebef) * htaft(j)) / dt hq(j) = ((timeaft - time) * hqbef(j) + (time - timebef) * hqaft(j)) / dt hw(j) = ((timeaft - time) * hwbef(j) + (time - timebef) * hwaft(j)) / dt IF(imp_fcg) THEN hu(j) = ((timeaft - time) * hubef(j) + (time - timebef) * huaft(j)) / dt hv(j) = ((timeaft - time) * hvbef(j) + (time - timebef) * hvaft(j)) / dt endif ! imp_fcg IF(Turb_fcg) THEN hThTurb(j) = ((timeaft - time) * hThTurbbef(j) & & + (time - timebef) * hThTurbaft(j)) / dt hqTurb(j) = ((timeaft - time) * hqTurbbef(j) & & + (time - timebef) * hqTurbaft(j)) / dt endif ! Turb_fcg enddo ts_subr = ((timeaft - time) * tsbef + (time - timebef) * tsaft) / dt endif ! temps.ge.pasmax print *, ' time,timebef,timeaft', time, timebef, timeaft print *, ' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft' do j = 1, klev print *, j, ht(j), htbef(j), htaft(j), & & hthturb(j), hthturbbef(j), hthturbaft(j) enddo print *, ' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft' do j = 1, klev print *, j, hq(j), hqbef(j), hqaft(j), & & hqturb(j), hqturbbef(j), hqturbaft(j) enddo !------------------------------------------------------------------- IF (Ts_fcg) Ts = Ts_subr return !----------------------------------------------------------------------- ! on sort les champs de "convergence" pour l instant initial 'in' ! ceci se passe au pas temps itap=0 de la physique !----------------------------------------------------------------------- entry get_uvd2(itap, dtime, file_forctl, file_fordat, & & ht, hq, hw, hu, hv, hThTurb, hqTurb, ts, & & imp_fcg, ts_fcg, Tp_fcg, Turb_fcg) PRINT*, 'le pas itap est:', itap !=================================================================== WRITE(*, '(a)') 'OPEN ' // file_forctl open(97, FILE = file_forctl, FORM = 'FORMATTED') !------------------ do i = 1, 1000 read(97, 1000, end = 999) string 1000 format (a4) IF (string == 'TDEF') go to 50 enddo 50 backspace(97) !------------------------------------------------------------------- ! *** on lit le pas de temps dans le fichier de donnees *** ! *** "forcing.ctl" et pasmax *** !------------------------------------------------------------------- read(97, 2000) aaa 2000 format (a80) PRINT*, 'aaa est', aaa aaa = spaces(aaa, 1) PRINT*, 'aaa', aaa CALL getsch(aaa, ' ', ' ', 5, atemps, nch) PRINT*, 'atemps est', atemps atemps = atemps(1:nch - 2) PRINT*, 'atemps', atemps read(atemps, *) imn dt = imn * 60 PRINT*, 'le pas de temps dt', dt CALL getsch(aaa, ' ', ' ', 2, apasmax, nch) apasmax = apasmax(1:nch) read(apasmax, *) ipa pasmax = ipa PRINT*, 'pasmax est', pasmax CLOSE(97) !------------------------------------------------------------------ ! *** on lit le pas de temps initial de la simulation *** !------------------------------------------------------------------ in = itap !c pasprev=in !c time0=dt*(pasprev-1) pasprev = in - 1 time0 = dt * pasprev close(98) WRITE(*, '(a)') 'OPEN ' // file_fordat open(99, FILE = file_fordat, FORM = 'UNFORMATTED', ACCESS = 'DIRECT', RECL = 8) icomp1 = nblvlm * 4 IF (ts_fcg) icomp1 = icomp1 + 1 IF (imp_fcg) icomp1 = icomp1 + nblvlm * 2 IF (Turb_fcg) icomp1 = icomp1 + nblvlm * 2 icompt = icomp1 * (in - 1) CALL rdgrads(99, icompt, nblvlm, z, ht_mes, hq_mes, hw_mes & &, hu_mes, hv_mes, hthturb_mes, hqturb_mes & &, ts_fcg, ts_subr, imp_fcg, Turb_fcg) print *, 'get_uvd : rdgrads ->' print *, tp_fcg ! following commented out because we have temperature already in ARM case ! (otherwise this is the potential temperature ) !------------------------------------------------------------------------ IF(Tp_fcg) THEN do i = 1, nblvlm ht_mes(i) = ht_mes(i) * (hplaym(i) / 1000.)**rkappa enddo endif ! Tp_fcg PRINT*, 'ht_mes ', (ht_mes(i), i = 1, nblvlm) PRINT*, 'hq_mes ', (hq_mes(i), i = 1, nblvlm) PRINT*, 'hw_mes ', (hw_mes(i), i = 1, nblvlm) IF(imp_fcg) THEN PRINT*, 'hu_mes ', (hu_mes(i), i = 1, nblvlm) PRINT*, 'hv_mes ', (hv_mes(i), i = 1, nblvlm) PRINT*, 't', ts_subr endif ! imp_fcg IF(Turb_fcg) THEN PRINT*, 'hThTurb_mes ', (hThTurb_mes(i), i = 1, nblvlm) PRINT*, 'hqTurb ', (hqTurb_mes(i), i = 1, nblvlm) endif ! Turb_fcg !---------------------------------------------------------------------- ! on a obtenu des champs initiaux sur les niveaux du meso_NH ! on interpole sur les niveaux du gcm(niveau pression bien sur!) !----------------------------------------------------------------------- do k = 1, klev IF (JM(k) == 0) THEN !FKC bug? ne faut il pas convertir tsol en tendance ???? !RT bug taken care of by removing the stuff htaft(k) = ht_mes(jm(k) + 1) hqaft(k) = hq_mes(jm(k) + 1) hwaft(k) = hw_mes(jm(k) + 1) IF(imp_fcg) THEN huaft(k) = hu_mes(jm(k) + 1) hvaft(k) = hv_mes(jm(k) + 1) endif ! imp_fcg IF(Turb_fcg) THEN hThTurbaft(k) = hThTurb_mes(jm(k) + 1) hqTurbaft(k) = hqTurb_mes(jm(k) + 1) endif ! Turb_fcg else ! JM(k) .EQ. 0 htaft(k) = coef1(k) * ht_mes(jm(k)) + coef2(k) * ht_mes(jm(k) + 1) hqaft(k) = coef1(k) * hq_mes(jm(k)) + coef2(k) * hq_mes(jm(k) + 1) hwaft(k) = coef1(k) * hw_mes(jm(k)) + coef2(k) * hw_mes(jm(k) + 1) IF(imp_fcg) THEN huaft(k) = coef1(k) * hu_mes(jm(k)) + coef2(k) * hu_mes(jm(k) + 1) hvaft(k) = coef1(k) * hv_mes(jm(k)) + coef2(k) * hv_mes(jm(k) + 1) endif ! imp_fcg IF(Turb_fcg) THEN hThTurbaft(k) = coef1(k) * hThTurb_mes(jm(k)) & & + coef2(k) * hThTurb_mes(jm(k) + 1) hqTurbaft(k) = coef1(k) * hqTurb_mes(jm(k)) & & + coef2(k) * hqTurb_mes(jm(k) + 1) endif ! Turb_fcg endif ! JM(k) .EQ. 0 enddo tsaft = ts_subr ! valeurs initiales des champs de convergence do k = 1, klev ht(k) = htaft(k) hq(k) = hqaft(k) hw(k) = hwaft(k) IF(imp_fcg) THEN hu(k) = huaft(k) hv(k) = hvaft(k) endif ! imp_fcg IF(Turb_fcg) THEN hThTurb(k) = hThTurbaft(k) hqTurb(k) = hqTurbaft(k) endif ! Turb_fcg enddo ts_subr = tsaft close(99) close(98) !------------------------------------------------------------------- IF (Ts_fcg) Ts = Ts_subr return 999 continue stop 'erreur lecture, file forcing.ctl' end SUBROUTINE advect_tvl(dtime, zt, zq, vu_f, vv_f, t_f, q_f & &, d_t_adv, d_q_adv) USE dimphy IMPLICIT NONE INCLUDE "dimensions.h" !cccc INCLUDE "dimphy.h" INTEGER k REAL dtime, fact, du, dv, cx, cy, alx, aly REAL zt(klev), zq(klev, 3) REAL vu_f(klev), vv_f(klev), t_f(klev), q_f(klev, 3) REAL d_t_adv(klev), d_q_adv(klev, 3) ! Velocity of moving cell data cx, cy /12., -2./ ! Dimensions of moving cell data alx, aly /100000., 150000./ do k = 1, klev du = abs(vu_f(k) - cx) / alx dv = abs(vv_f(k) - cy) / aly fact = dtime * (du + dv - du * dv * dtime) d_t_adv(k) = fact * (t_f(k) - zt(k)) d_q_adv(k, 1) = fact * (q_f(k, 1) - zq(k, 1)) d_q_adv(k, 2) = fact * (q_f(k, 2) - zq(k, 2)) d_q_adv(k, 3) = fact * (q_f(k, 3) - zq(k, 3)) enddo RETURN end SUBROUTINE copie(klevgcm, playgcm, psolgcm, file_forctl) IMPLICIT NONE !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER klev !nombre de niveau de pression du GCM REAL play(100) !pression en Pa au milieu de chaque couche GCM INTEGER JM(100) REAL coef1(100) !coefficient d interpolation REAL coef2(100) !coefficient d interpolation INTEGER nblvlm !nombre de niveau de pression du mesoNH REAL playm(100) !pression en Pa au milieu de chaque couche Meso-NH REAL hplaym(100)!pression en hecto-Pa des milieux de couche Meso-NH COMMON/com1_phys_gcss/play, coef1, coef2, JM, klev COMMON/com2_phys_gcss/playm, hplaym, nblvlm INTEGER k, klevgcm REAL playgcm(klevgcm) ! pression en milieu de couche du gcm REAL psolgcm CHARACTER*80 file_forctl klev = klevgcm !--------------------------------------------------------------------- ! pression au milieu des couches du gcm dans la physiq ! (SB: remplace le CALL conv_lipress_gcm(playgcm) ) !--------------------------------------------------------------------- do k = 1, klev play(k) = playgcm(k) PRINT*, 'la pression gcm est:', play(k) enddo !---------------------------------------------------------------------- ! lecture du descripteur des donnees Meso-NH (forcing.ctl): ! -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH ! (on remplit le COMMON com2_phys_gcss) !---------------------------------------------------------------------- CALL mesolupbis(file_forctl) PRINT*, 'la valeur de nblvlm est:', nblvlm !---------------------------------------------------------------------- ! etude de la correspondance entre les niveaux meso.NH et GCM; ! calcul des coefficients d interpolation coef1 et coef2 ! (on remplit le COMMON com1_phys_gcss) !---------------------------------------------------------------------- CALL corresbis(psolgcm) !--------------------------------------------------------- ! TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss: !--------------------------------------------------------- WRITE(*, *) ' ' WRITE(*, *) 'TESTS com1_phys_gcss et com2_phys_gcss dans copie.F' WRITE(*, *) '--------------------------------------' WRITE(*, *) 'GCM: nb niveaux:', klev, ' et pression, coeffs:' do k = 1, klev WRITE(*, *) play(k), coef1(k), coef2(k) enddo WRITE(*, *) 'MESO-NH: nb niveaux:', nblvlm, ' et pression:' do k = 1, nblvlm WRITE(*, *) playm(k), hplaym(k) enddo WRITE(*, *) ' ' END SUBROUTINE mesolupbis(file_forctl) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Lecture descripteur des donnees MESO-NH (forcing.ctl): ! ------------------------------------------------------- ! Cette subroutine lit dans le fichier de controle "essai.ctl" ! et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs ! des pressions en milieu de couche du Meso-NH (en Pa puis en hPa). !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER nblvlm !nombre de niveau de pression du mesoNH REAL playm(100) !pression en Pa milieu de chaque couche Meso-NH REAL hplaym(100) !pression en hPa des milieux de couche Meso-NH COMMON/com2_phys_gcss/playm, hplaym, nblvlm INTEGER i, lu, mlz, mlzh CHARACTER*80 file_forctl CHARACTER*4 a CHARACTER*80 aaa, anblvl INTEGER nch lu = 9 open(lu, file = file_forctl, form = 'formatted') do i = 1, 1000 read(lu, 1000, end = 999) a IF (a == 'ZDEF') go to 100 enddo 100 backspace(lu) PRINT*, ' DESCRIPTION DES 2 MODELES : ' PRINT*, ' ' read(lu, 2000) aaa 2000 format (a80) aaa = spaces(aaa, 1) CALL getsch(aaa, ' ', ' ', 2, anblvl, nch) read(anblvl, *) nblvlm PRINT*, 'nbre de niveaux de pression Meso-NH :', nblvlm PRINT*, ' ' PRINT*, 'pression en Pa de chaque couche du meso-NH :' read(lu, *) (playm(mlz), mlz = 1, nblvlm) ! Si la pression est en HPa, la multiplier par 100 IF (playm(1) < 10000.) THEN do mlz = 1, nblvlm playm(mlz) = playm(mlz) * 100. enddo endif PRINT*, (playm(mlz), mlz = 1, nblvlm) 1000 format (a4) PRINT*, ' ' do mlzh = 1, nblvlm hplaym(mlzh) = playm(mlzh) / 100. enddo PRINT*, 'pression en hPa de chaque couche du meso-NH: ' PRINT*, (hplaym(mlzh), mlzh = 1, nblvlm) close (lu) return 999 stop 'erreur lecture des niveaux pression des donnees' end SUBROUTINE rdgrads(itape, icount, nl, z, ht, hq, hw, hu, hv, hthtur, hqtur, & & ts_fcg, ts, imp_fcg, Turb_fcg) IMPLICIT NONE INTEGER itape, icount, icomp, nl REAL z(nl), ht(nl), hq(nl), hw(nl), hu(nl), hv(nl) REAL hthtur(nl), hqtur(nl) REAL ts INTEGER k LOGICAL imp_fcg, ts_fcg, Turb_fcg icomp = icount do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)z(k) print *, 'icomp,k,z(k) ', icomp, k, z(k) enddo do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)hT(k) PRINT*, hT(k), k enddo do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)hQ(k) enddo IF(turb_fcg) THEN do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)hThTur(k) enddo do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)hqTur(k) enddo endif print *, ' apres lecture hthtur, hqtur' IF(imp_fcg) THEN do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)hu(k) enddo do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)hv(k) enddo endif do k = 1, nl icomp = icomp + 1 read(itape, rec = icomp)hw(k) enddo IF(ts_fcg) THEN icomp = icomp + 1 read(itape, rec = icomp)ts endif print *, ' rdgrads ->' RETURN END SUBROUTINE corresbis(psol) IMPLICIT NONE !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Cette subroutine calcule et affiche les valeurs des coefficients ! d interpolation qui serviront dans la formule d interpolation elle- ! meme. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER klev !nombre de niveau de pression du GCM REAL play(100) !pression en Pa au milieu de chaque couche GCM INTEGER JM(100) REAL coef1(100) !coefficient d interpolation REAL coef2(100) !coefficient d interpolation INTEGER nblvlm !nombre de niveau de pression du mesoNH REAL playm(100) !pression en Pa milieu de chaque couche Meso-NH REAL hplaym(100)!pression en hPa des milieux de couche Meso-NH COMMON/com1_phys_gcss/play, coef1, coef2, JM, klev COMMON/com2_phys_gcss/playm, hplaym, nblvlm REAL psol REAL val INTEGER k, mlz do k = 1, klev val = play(k) IF (val > playm(1)) THEN mlz = 0 JM(1) = mlz coef1(1) = (playm(mlz + 1) - val) / (playm(mlz + 1) - psol) coef2(1) = (val - psol) / (playm(mlz + 1) - psol) ELSE IF (val > playm(nblvlm)) THEN do mlz = 1, nblvlm IF (val <= playm(mlz).AND. val > playm(mlz + 1))THEN JM(k) = mlz coef1(k) = (playm(mlz + 1) - val) / (playm(mlz + 1) - playm(mlz)) coef2(k) = (val - playm(mlz)) / (playm(mlz + 1) - playm(mlz)) endif enddo else JM(k) = nblvlm - 1 coef1(k) = 0. coef2(k) = 0. endif enddo !c if (play(klev) .le. playm(nblvlm)) THEN !c mlz=nblvlm-1 !c JM(klev)=mlz !c coef1(klev)=(playm(mlz+1)-val) !c * /(playm(mlz+1)-playm(mlz)) !c coef2(klev)=(val-playm(mlz)) !c * /(playm(mlz+1)-playm(mlz)) !c endif PRINT*, ' ' PRINT*, ' INTERPOLATION : ' PRINT*, ' ' PRINT*, 'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:' PRINT*, (JM(k), k = 1, klev) PRINT*, 'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:' PRINT*, (JM(k), k = 1, klev) PRINT*, ' ' PRINT*, 'vals du premier coef d"interpolation pour les 9 niveaux: ' PRINT*, (coef1(k), k = 1, klev) PRINT*, ' ' PRINT*, 'valeurs du deuxieme coef d"interpolation pour les 9 niveaux:' PRINT*, (coef2(k), k = 1, klev) RETURN END SUBROUTINE GETSCH(STR, DEL, TRM, NTH, SST, NCH) !*************************************************************** !* * !* * !* GETSCH * !* * !* * !* modified by : * !*************************************************************** !* Return in SST the character string found between the NTH-1 and NTH !* occurence of the delimiter 'DEL' but before the terminator 'TRM' in !* the input string 'STR'. If TRM=DEL then STR is considered unlimited. !* NCH=Length of the string returned in SST or =-1 if NTH is <1 or if !* NTH is greater than the number of delimiters in STR. IMPLICIT INTEGER (A-Z) CHARACTER STR*(*), DEL*1, TRM*1, SST*(*) NCH = -1 SST = ' ' IF(NTH>0) THEN IF(TRM==DEL) THEN LENGTH = LEN(STR) ELSE LENGTH = INDEX(STR, TRM) - 1 IF(LENGTH<0) LENGTH = LEN(STR) ENDIF !* Find beginning and end of the NTH DEL-limited substring in STR END = -1 DO N = 1, NTH IF(END==LENGTH) RETURN BEG = END + 2 END = BEG + INDEX(STR(BEG:LENGTH), DEL) - 2 IF(END==BEG - 2) END = LENGTH !* PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END end DO NCH = END - BEG + 1 IF(NCH>0) SST = STR(BEG:END) ENDIF END CHARACTER*(80) FUNCTION SPACES(STR, NSPACE) ! CERN PROGLIB# M433 SPACES .VERSION KERNFOR 4.14 860211 ! ORIG. 6/05/86 M.GOOSSENS/DD !- The function value SPACES returns the character string STR with !- leading blanks removed and each occurence of one or more blanks !- replaced by NSPACE blanks inside the string STR CHARACTER*(80) STR INTEGER nspace, IBLANK, ISPACE, INONBL, LENSPA LENSPA = LEN(SPACES) SPACES = ' ' IF (NSPACE<0) NSPACE = 0 IBLANK = 1 ISPACE = 1 100 INONBL = INDEXC(STR(IBLANK:), ' ') IF (INONBL==0) THEN SPACES(ISPACE:) = STR(IBLANK:) RETURN ENDIF INONBL = INONBL + IBLANK - 1 IBLANK = INDEX(STR(INONBL:), ' ') IF (IBLANK==0) THEN SPACES(ISPACE:) = STR(INONBL:) RETURN ENDIF IBLANK = IBLANK + INONBL - 1 SPACES(ISPACE:) = STR(INONBL:IBLANK - 1) ISPACE = ISPACE + IBLANK - INONBL + NSPACE IF (ISPACE<=LENSPA) GO TO 100 END INTEGER FUNCTION INDEXC(STR, SSTR) ! CERN PROGLIB# M433 INDEXC .VERSION KERNFOR 4.14 860211 ! ORIG. 26/03/86 M.GOOSSENS/DD !- Find the leftmost position where substring SSTR does not match !- string STR scanning forward CHARACTER*(*) STR, SSTR INTEGER I, LENS, LENSS LENS = LEN(STR) LENSS = LEN(SSTR) DO I = 1, LENS - LENSS + 1 IF (STR(I:I + LENSS - 1)/=SSTR) THEN INDEXC = I RETURN ENDIF END DO INDEXC = 0 END END MODULE lmdz_old_1dconv