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) c implicit none cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de c pouvoir calculer la convergence et le cisaillement dans la physiq ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc #include "YOMCST.h" 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,ii,ll,in REAL tsol,qsol CHARACTER*80 file_forctl,file_fordat COMMON/com1_phys_gcss/klev,play,JM,coef1,coef2 COMMON/com2_phys_gcss/nblvlm,playm,hplaym c====================================================================== c methode: on va chercher les donnees du mesoNH de meteo france, on y c a acces a tout pas detemps grace a la routine rdgrads qui c est une boucle lisant dans ces fichiers. c Puis on interpole ces donnes sur les 11 niveaux du gcm et c et sur les pas de temps de ce meme gcm c---------------------------------------------------------------------- c input: c pasmax :nombre de pas de temps maximum du mesoNH c dt :pas de temps du meso_NH (en secondes) c---------------------------------------------------------------------- integer pasmax,dt save pasmax,dt c---------------------------------------------------------------------- c arguments: c itap :compteur de la physique(le nombre de ces pas est c fixe dans la subroutine calcul_ini_gcm de interpo c -lation c dtime :pas detemps du gcm (en secondes) c ht :convergence horizontale de temperature(K/s) c hq : " " d humidite (kg/kg/s) c hw :vitesse verticale moyenne (m/s**2) c hu :convergence horizontale d impulsion le long de x c (kg/(m^2 s^2) c hv : idem le long de y. c Ts : Temperature de surface (K) c imp_fcg: var. logical .eq. T si forcage en impulsion c ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier c Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle c Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier c---------------------------------------------------------------------- integer itap real dtime real ht(100) real hq(100) real hu(100) real hv(100) real hw(100) real hthturb(100) real hqturb(100) real Ts, Ts_subr logical imp_fcg logical ts_fcg logical Tp_fcg logical Turb_fcg c---------------------------------------------------------------------- c Variables internes de get_uvd (note : l interpolation temporelle c est faite entre les pas de temps before et after, sur les variables c definies sur la grille du SCM; on atteint exactement les valeurs Meso c aux milieux des pas de temps Meso) c time0 :date initiale en secondes c time :temps associe a chaque pas du SCM c pas :numero du pas du meso_NH (on lit en pas : le premier pas c des donnees est duplique) c pasprev :numero du pas de lecture precedent c htaft :advection horizontale de temp. au pas de temps after c hqaft : " " d humidite " c hwaft :vitesse verticalle moyenne au pas de temps after c huaft,hvaft :advection horizontale d impulsion au pas de temps after c tsaft : surface temperature 'after time step' c htbef :idem htaft, mais pour le pas de temps before c hqbef :voir hqaft c hwbef :voir hwaft c hubef,hvbef : idem huaft,hvaft, mais pour before c tsbef : surface temperature 'before time step' c---------------------------------------------------------------------- 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 c real timeaft,timebef save timeaft,timebef integer temps character*4 string c---------------------------------------------------------------------- c variables arguments de la subroutine rdgrads c--------------------------------------------------------------------- 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 c c--------------------------------------------------------------------- c variable argument de la subroutine copie c--------------------------------------------------------------------- c SB real pplay(100) !pression en milieu de couche du gcm c SB !argument de la physique c--------------------------------------------------------------------- c variables destinees a la lecture du pas de temps du fichier de donnees c--------------------------------------------------------------------- character*80 aaa,atemps,spaces,apasmax integer nch,imn,ipa c--------------------------------------------------------------------- c procedures appelees external rdgrads !lire en iterant dans forcing.dat c--------------------------------------------------------------------- print*,'le pas itap est:',itap c*** on determine le pas du meso_NH correspondant au nouvel itap *** c*** pour aller chercher les champs dans rdgrads *** c time=time0+itap*dtime cc temps=int(time/dt+1) cc pas=min(temps,pasmax) temps = 1 + int((dt + 2*time)/(2*dt)) pas=min(temps,pasmax-1) print*,'le pas Meso est:',pas c c c=================================================================== c c*** on remplit les champs before avec les champs after du pas *** c*** precedent en format gcm *** if(pas.gt.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 c*** on va chercher les nouveaux champs after dans toga.dat *** c*** 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) c if(Tp_fcg) then c (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 c 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 c*** on interpole les champs meso_NH sur les niveaux de pression*** c*** gcm . on obtient le nouveau champ after *** do k=1,klev if (JM(k) .eq. 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 c*** si on atteint le pas max des donnees experimentales ,on *** c*** on conserve les derniers champs calcules *** if(temps.ge.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 c*** on interpole sur les pas de temps de 10mn du gcm a partir *** c** 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 c 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 c c------------------------------------------------------------------- c IF (Ts_fcg) Ts = Ts_subr return c c----------------------------------------------------------------------- c on sort les champs de "convergence" pour l instant initial 'in' c ceci se passe au pas temps itap=0 de la physique c----------------------------------------------------------------------- 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 c c=================================================================== c write(*,'(a)') 'OPEN '//file_forctl open(97,FILE=file_forctl,FORM='FORMATTED') c c------------------ do i=1,1000 read(97,1000,end=999) string 1000 format (a4) if (string .eq. 'TDEF') go to 50 enddo 50 backspace(97) c------------------------------------------------------------------- c *** on lit le pas de temps dans le fichier de donnees *** c *** "forcing.ctl" et pasmax *** c------------------------------------------------------------------- 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) c------------------------------------------------------------------ c *** on lit le pas de temps initial de la simulation *** c------------------------------------------------------------------ in=itap cc pasprev=in cc time0=dt*(pasprev-1) pasprev=in-1 time0=dt*pasprev C close(98) c 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 c c following commented out because we have temperature already in ARM case c (otherwise this is the potential temperature ) c------------------------------------------------------------------------ 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 c---------------------------------------------------------------------- c on a obtenu des champs initiaux sur les niveaux du meso_NH c on interpole sur les niveaux du gcm(niveau pression bien sur!) c----------------------------------------------------------------------- do k=1,klev if (JM(k) .eq. 0) then cFKC bug? ne faut il pas convertir tsol en tendance ???? cRT 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 c 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) c c------------------------------------------------------------------- c c 100 IF (Ts_fcg) Ts = Ts_subr return c 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" ccccc#include "dimphy.h" integer k real dtime, fact, du, dv, cx, cy, alx, aly real zt(klev), zq(klev,3) : , vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3) real d_t_adv(klev), d_q_adv(klev,3) c Velocity of moving cell data cx,cy /12., -2./ c Dimensions of moving cell data alx,aly /100 000.,150 000./ 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 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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/klev,play,JM,coef1,coef2 COMMON/com2_phys_gcss/nblvlm,playm,hplaym integer i,k,klevgcm real playgcm(klevgcm) ! pression en milieu de couche du gcm real psolgcm character*80 file_forctl klev = klevgcm c--------------------------------------------------------------------- c pression au milieu des couches du gcm dans la physiq c (SB: remplace le call conv_lipress_gcm(playgcm) ) c--------------------------------------------------------------------- do k = 1, klev play(k) = playgcm(k) print*,'la pression gcm est:',play(k) enddo c---------------------------------------------------------------------- c lecture du descripteur des donnees Meso-NH (forcing.ctl): c -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH c (on remplit le COMMON com2_phys_gcss) c---------------------------------------------------------------------- call mesolupbis(file_forctl) print*,'la valeur de nblvlm est:',nblvlm c---------------------------------------------------------------------- c etude de la correspondance entre les niveaux meso.NH et GCM; c calcul des coefficients d interpolation coef1 et coef2 c (on remplit le COMMON com1_phys_gcss) c---------------------------------------------------------------------- call corresbis(psolgcm) c--------------------------------------------------------- c TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss: c--------------------------------------------------------- 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) implicit none c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Lecture descripteur des donnees MESO-NH (forcing.ctl): c ------------------------------------------------------- c c Cette subroutine lit dans le fichier de controle "essai.ctl" c et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs c des pressions en milieu de couche du Meso-NH (en Pa puis en hPa). cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 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/nblvlm,playm,hplaym INTEGER i,lu,k,mlz,mlzh,j character*80 file_forctl character*4 a character*80 aaa,anblvl,spaces integer nch lu=9 open(lu,file=file_forctl,form='formatted') c do i=1,1000 read(lu,1000,end=999) a if (a .eq. 'ZDEF') go to 100 enddo c 100 backspace(lu) print*,' DESCRIPTION DES 2 MODELES : ' print*,' ' c read(lu,2000) aaa 2000 format (a80) aaa=spaces(aaa,1) call getsch(aaa,' ',' ',2,anblvl,nch) read(anblvl,*) nblvlm c print*,'nbre de niveaux de pression Meso-NH :',nblvlm print*,' ' print*,'pression en Pa de chaque couche du meso-NH :' c read(lu,*) (playm(mlz),mlz=1,nblvlm) c Si la pression est en HPa, la multiplier par 100 if (playm(1) .lt. 10000.) then do mlz = 1,nblvlm playm(mlz) = playm(mlz)*100. enddo endif print*,(playm(mlz),mlz=1,nblvlm) c 1000 format (a4) 1001 format(5x,i2) c print*,' ' do mlzh=1,nblvlm hplaym(mlzh)=playm(mlzh)/100. enddo c print*,'pression en hPa de chaque couche du meso-NH: ' print*,(hplaym(mlzh),mlzh=1,nblvlm) c close (lu) return c 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 c INTEGER i, k c LOGICAL imp_fcg,ts_fcg,Turb_fcg c icomp = icount c c 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 c 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' c 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 c do k=1,nl icomp=icomp+1 read(itape,rec=icomp)hw(k) enddo c if(ts_fcg) then icomp=icomp+1 read(itape,rec=icomp)ts endif c print *,' rdgrads ->' RETURN END SUBROUTINE corresbis(psol) implicit none ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Cette subroutine calcule et affiche les valeurs des coefficients c d interpolation qui serviront dans la formule d interpolation elle- c meme. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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/klev,play,JM,coef1,coef2 COMMON/com2_phys_gcss/nblvlm,playm,hplaym REAL psol REAL val INTEGER k, mlz, mlzh do k=1,klev val=play(k) if (val .gt. 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 .gt. playm(nblvlm)) then do mlz=1,nblvlm if ( val .le. playm(mlz) * .and. val .gt. 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 cc if (play(klev) .le. playm(nblvlm)) then cc mlz=nblvlm-1 cc JM(klev)=mlz cc coef1(klev)=(playm(mlz+1)-val) cc * /(playm(mlz+1)-playm(mlz)) cc coef2(klev)=(val-playm(mlz)) cc * /(playm(mlz+1)-playm(mlz)) cc endif c 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*,'valeurs 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 niveau *x: ' print*,(coef2(k),k=1,klev) c return end SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH) C*************************************************************** C* * C* * C* GETSCH * C* * C* * C* modified by : * C*************************************************************** C* Return in SST the character string found between the NTH-1 and NTH C* occurence of the delimiter 'DEL' but before the terminator 'TRM' in C* the input string 'STR'. If TRM=DEL then STR is considered unlimited. C* NCH=Length of the string returned in SST or =-1 if NTH is <1 or if C* 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.GT.0) THEN IF(TRM.EQ.DEL) THEN LENGTH=LEN(STR) ELSE LENGTH=INDEX(STR,TRM)-1 IF(LENGTH.LT.0) LENGTH=LEN(STR) ENDIF C* Find beginning and end of the NTH DEL-limited substring in STR END=-1 DO 1,N=1,NTH IF(END.EQ.LENGTH) RETURN BEG=END+2 END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2 IF(END.EQ.BEG-2) END=LENGTH C* PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END 1 CONTINUE NCH=END-BEG+1 IF(NCH.GT.0) SST=STR(BEG:END) ENDIF END CHARACTER*(*) FUNCTION SPACES(STR,NSPACE) C C CERN PROGLIB# M433 SPACES .VERSION KERNFOR 4.14 860211 C ORIG. 6/05/86 M.GOOSSENS/DD C C- The function value SPACES returns the character string STR with C- leading blanks removed and each occurence of one or more blanks C- replaced by NSPACE blanks inside the string STR C CHARACTER*(*) STR C LENSPA = LEN(SPACES) SPACES = ' ' IF (NSPACE.LT.0) NSPACE = 0 IBLANK = 1 ISPACE = 1 100 INONBL = INDEXC(STR(IBLANK:),' ') IF (INONBL.EQ.0) THEN SPACES(ISPACE:) = STR(IBLANK:) GO TO 999 ENDIF INONBL = INONBL + IBLANK - 1 IBLANK = INDEX(STR(INONBL:),' ') IF (IBLANK.EQ.0) THEN SPACES(ISPACE:) = STR(INONBL:) GO TO 999 ENDIF IBLANK = IBLANK + INONBL - 1 SPACES(ISPACE:) = STR(INONBL:IBLANK-1) ISPACE = ISPACE + IBLANK - INONBL + NSPACE IF (ISPACE.LE.LENSPA) GO TO 100 999 END FUNCTION INDEXC(STR,SSTR) C C CERN PROGLIB# M433 INDEXC .VERSION KERNFOR 4.14 860211 C ORIG. 26/03/86 M.GOOSSENS/DD C C- Find the leftmost position where substring SSTR does not match C- string STR scanning forward C CHARACTER*(*) STR,SSTR C LENS = LEN(STR) LENSS = LEN(SSTR) C DO 10 I=1,LENS-LENSS+1 IF (STR(I:I+LENSS-1).NE.SSTR) THEN INDEXC = I GO TO 999 ENDIF 10 CONTINUE INDEXC = 0 C 999 END