! ! $Id: conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec $ ! c c SUBROUTINE conf_unicol( tapedef ) c #ifdef CPP_IOIPSL use IOIPSL #else ! if not using IOIPSL, we still need to use (a local version of) getin use ioipsl_getincom #endif IMPLICIT NONE c----------------------------------------------------------------------- c Auteurs : A. Lahellec . c c Arguments : c c tapedef : INTEGER tapedef c c Declarations : c -------------- #include "compar1d.h" #include "flux_arp.h" #include "tsoilnudge.h" #include "fcg_gcssold.h" #include "fcg_racmo.h" #include "iniprint.h" c c c local: c ------ c CHARACTER ch1*72,ch2*72,ch3*72,ch4*12 c c ------------------------------------------------------------------- c c ......... Initilisation parametres du lmdz1D .......... c c--------------------------------------------------------------------- c initialisations: c ---------------- !Config Key = lunout !Config Desc = unite de fichier pour les impressions !Config Def = 6 !Config Help = unite de fichier pour les impressions !Config (defaut sortie standard = 6) lunout=6 ! CALL getin('lunout', lunout) IF (lunout /= 5 .and. lunout /= 6) THEN OPEN(lunout,FILE='lmdz.out') ENDIF !Config Key = prt_level !Config Desc = niveau d'impressions de débogage !Config Def = 0 !Config Help = Niveau d'impression pour le débogage !Config (0 = minimum d'impression) c prt_level = 0 c CALL getin('prt_level',prt_level) c----------------------------------------------------------------------- c Parametres de controle du run: c----------------------------------------------------------------------- !Config Key = restart !Config Desc = on repart des startphy et start1dyn !Config Def = false !Config Help = les fichiers restart doivent etre renomme en start restart = .FALSE. CALL getin('restart',restart) !Config Key = forcing_type !Config Desc = defines the way the SCM is forced: !Config Def = 0 !!Config Help = 0 ==> forcing_les = .true. ! initial profiles from file prof.inp.001 ! no forcing by LS convergence ; ! surface temperature imposed ; ! radiative cooling may be imposed (iflag_radia=0 in physiq.def) ! = 1 ==> forcing_radconv = .true. ! idem forcing_type = 0, but the imposed radiative cooling ! is set to 0 (hence, if iflag_radia=0 in physiq.def, ! then there is no radiative cooling at all) ! = 2 ==> forcing_toga = .true. ! initial profiles from TOGA-COARE IFA files ! LS convergence and SST imposed from TOGA-COARE IFA files ! = 3 ==> forcing_GCM2SCM = .true. ! initial profiles from the GCM output ! LS convergence imposed from the GCM output ! = 4 ==> forcing_twpi = .true. ! initial profiles from TWPICE nc files ! LS convergence and SST imposed from TWPICE nc files ! = 5 ==> forcing_rico = .true. ! initial profiles from RICO idealized ! LS convergence imposed from RICO (cst) ! = 6 ==> forcing_amma = .true. ! = 40 ==> forcing_GCSSold = .true. ! initial profile from GCSS file ! LS convergence imposed from GCSS file ! = 59 ==> forcing_sandu = .true. ! initial profiles from sanduref file: see prof.inp.001 ! SST varying with time and divergence constante: see ifa_sanduref.txt file ! Radiation has to be computed interactively ! = 60 ==> forcing_astex = .true. ! initial profiles from file: see prof.inp.001 ! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file ! Radiation has to be computed interactively ! = 61 ==> forcing_armcu = .true. ! initial profiles from file: see prof.inp.001 ! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt ! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt ! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s ! Radiation to be switched off ! forcing_type = 0 CALL getin('forcing_type',forcing_type) imp_fcg_gcssold = .false. ts_fcg_gcssold = .false. Tp_fcg_gcssold = .false. Tp_ini_gcssold = .false. xTurb_fcg_gcssold = .false. IF (forcing_type .eq.40) THEN CALL getin('imp_fcg',imp_fcg_gcssold) CALL getin('ts_fcg',ts_fcg_gcssold) CALL getin('tp_fcg',Tp_fcg_gcssold) CALL getin('tp_ini',Tp_ini_gcssold) CALL getin('turb_fcg',xTurb_fcg_gcssold) ENDIF !Config Key = ok_flux_surf !Config Desc = forcage ou non par les flux de surface !Config Def = false !Config Help = forcage ou non par les flux de surface ok_flux_surf = .FALSE. CALL getin('ok_flux_surf',ok_flux_surf) !Config Key = ok_old_disvert !Config Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h) !Config Def = false !Config Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h) ok_old_disvert = .FALSE. CALL getin('ok_old_disvert',ok_old_disvert) !Config Key = time_ini !Config Desc = meaningless in this case !Config Def = 0. !Config Help = tsurf = 0. CALL getin('time_ini',time_ini) !Config Key = rlat et rlon !Config Desc = latitude et longitude !Config Def = 0.0 0.0 !Config Help = fixe la position de la colonne xlat = 0. xlon = 0. CALL getin('rlat',xlat) CALL getin('rlon',xlon) !Config Key = airephy !Config Desc = Grid cell area !Config Def = 1.e11 !Config Help = airefi = 1.e11 CALL getin('airephy',airefi) !Config Key = nat_surf !Config Desc = surface type !Config Def = 0 (ocean) !Config Help = 0=ocean,1=land,2=glacier,3=banquise nat_surf = 0. CALL getin('nat_surf',nat_surf) !Config Key = tsurf !Config Desc = surface temperature !Config Def = 290. !Config Help = not used if type_ts_forcing=1 in lmdz1d.F tsurf = 290. CALL getin('tsurf',tsurf) !Config Key = psurf !Config Desc = surface pressure !Config Def = 102400. !Config Help = psurf = 102400. CALL getin('psurf',psurf) !Config Key = zsurf !Config Desc = surface altitude !Config Def = 0. !Config Help = zsurf = 0. CALL getin('zsurf',zsurf) !Config Key = rugos !Config Desc = coefficient de frottement !Config Def = 0.0001 !Config Help = calcul du Cdrag rugos = 0.0001 CALL getin('rugos',rugos) !Config Key = wtsurf et wqsurf !Config Desc = ??? !Config Def = 0.0 0.0 !Config Help = wtsurf = 0.0 wqsurf = 0.0 CALL getin('wtsurf',wtsurf) CALL getin('wqsurf',wqsurf) !Config Key = albedo !Config Desc = albedo !Config Def = 0.09 !Config Help = albedo = 0.09 CALL getin('albedo',albedo) !Config Key = agesno !Config Desc = age de la neige !Config Def = 30.0 !Config Help = xagesno = 30.0 CALL getin('agesno',xagesno) !Config Key = restart_runoff !Config Desc = age de la neige !Config Def = 30.0 !Config Help = restart_runoff = 0.0 CALL getin('restart_runoff',restart_runoff) !Config Key = qsolinp !Config Desc = initial bucket water content (kg/m2) when land (5std) !Config Def = 30.0 !Config Help = qsolinp = 1. CALL getin('qsolinp',qsolinp) !Config Key = zpicinp !Config Desc = denivellation orographie !Config Def = 300. !Config Help = input brise zpicinp = 300. CALL getin('zpicinp',zpicinp) !Config key = nudge_tsoil !Config Desc = activation of soil temperature nudging !Config Def = .FALSE. !Config Help = ... nudge_tsoil=.FALSE. CALL getin('nudge_tsoil',nudge_tsoil) !Config key = isoil_nudge !Config Desc = level number where soil temperature is nudged !Config Def = 3 !Config Help = ... isoil_nudge=3 CALL getin('isoil_nudge',isoil_nudge) !Config key = Tsoil_nudge !Config Desc = target temperature for tsoil(isoil_nudge) !Config Def = 300. !Config Help = ... Tsoil_nudge=300. CALL getin('Tsoil_nudge',Tsoil_nudge) !Config key = tau_soil_nudge !Config Desc = nudging relaxation time for tsoil !Config Def = 3600. !Config Help = ... tau_soil_nudge=3600. CALL getin('tau_soil_nudge',tau_soil_nudge) write(lunout,*)' +++++++++++++++++++++++++++++++++++++++' write(lunout,*)' Configuration des parametres du gcm1D: ' write(lunout,*)' +++++++++++++++++++++++++++++++++++++++' write(lunout,*)' restart = ', restart write(lunout,*)' forcing_type = ', forcing_type write(lunout,*)' time_ini = ', time_ini write(lunout,*)' rlat = ', xlat write(lunout,*)' rlon = ', xlon write(lunout,*)' airephy = ', airefi write(lunout,*)' nat_surf = ', nat_surf write(lunout,*)' tsurf = ', tsurf write(lunout,*)' psurf = ', psurf write(lunout,*)' zsurf = ', zsurf write(lunout,*)' rugos = ', rugos write(lunout,*)' wtsurf = ', wtsurf write(lunout,*)' wqsurf = ', wqsurf write(lunout,*)' albedo = ', albedo write(lunout,*)' xagesno = ', xagesno write(lunout,*)' restart_runoff = ', restart_runoff write(lunout,*)' qsolinp = ', qsolinp write(lunout,*)' zpicinp = ', zpicinp write(lunout,*)' nudge_tsoil = ', nudge_tsoil write(lunout,*)' isoil_nudge = ', isoil_nudge write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge IF (forcing_type .eq.40) THEN write(lunout,*) '--- Forcing type GCSS Old --- with:' write(lunout,*)'imp_fcg',imp_fcg_gcssold write(lunout,*)'ts_fcg',ts_fcg_gcssold write(lunout,*)'tp_fcg',Tp_fcg_gcssold write(lunout,*)'tp_ini',Tp_ini_gcssold write(lunout,*)'xturb_fcg',xTurb_fcg_gcssold ENDIF write(lunout,*)' +++++++++++++++++++++++++++++++++++++++' write(lunout,*) c RETURN END ! ! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$ ! c SUBROUTINE dyn1deta0(fichnom,plev,play,phi,phis,presnivs, & ucov,vcov,temp,q,omega2) USE dimphy USE mod_grid_phy_lmdz USE mod_phys_lmdz_para USE iophy USE phys_state_var_mod USE iostart USE write_field_phy USE infotrac use control_mod IMPLICIT NONE c======================================================= c Ecriture du fichier de redemarrage sous format NetCDF c======================================================= c Declarations: c ------------- #include "dimensions.h" #include "comconst.h" #include "temps.h" !!#include "control.h" #include "logic.h" #include "netcdf.inc" c Arguments: c ---------- CHARACTER*(*) fichnom cAl1 plev tronque pour .nc mais plev(klev+1):=0 real :: plev(klon,klev),play (klon,klev),phi(klon,klev) real :: presnivs(klon,klev) real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev) real :: q(klon,klev,nqtot),omega2(klon,klev) c real :: ug(klev),vg(klev),fcoriolis real :: phis(klon) c Variables locales pour NetCDF: c ------------------------------ INTEGER nid, nvarid INTEGER idim_s INTEGER ierr, ierr_file INTEGER iq INTEGER length PARAMETER (length = 100) REAL tab_cntrl(length) ! tableau des parametres du run character*4 nmq(nqtot) character*12 modname character*80 abort_message LOGICAL found c INTEGER nb modname = 'dyn1deta0 : ' nmq(1)="vap" nmq(2)="cond" do iq=3,nqtot write(nmq(iq),'("tra",i1)') iq-2 enddo print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot CALL open_startphy(fichnom) print*,'after open startphy ',fichnom,nmq c c Lecture des parametres de controle: c CALL get_var("controle",tab_cntrl) im = tab_cntrl(1) jm = tab_cntrl(2) lllm = tab_cntrl(3) day_ref = tab_cntrl(4) annee_ref = tab_cntrl(5) c rad = tab_cntrl(6) c omeg = tab_cntrl(7) c g = tab_cntrl(8) c cpp = tab_cntrl(9) c kappa = tab_cntrl(10) c daysec = tab_cntrl(11) c dtvr = tab_cntrl(12) c etot0 = tab_cntrl(13) c ptot0 = tab_cntrl(14) c ztot0 = tab_cntrl(15) c stot0 = tab_cntrl(16) c ang0 = tab_cntrl(17) c pa = tab_cntrl(18) c preff = tab_cntrl(19) c c clon = tab_cntrl(20) c clat = tab_cntrl(21) c grossismx = tab_cntrl(22) c grossismy = tab_cntrl(23) c IF ( tab_cntrl(24).EQ.1. ) THEN fxyhypb = . TRUE . c dzoomx = tab_cntrl(25) c dzoomy = tab_cntrl(26) c taux = tab_cntrl(28) c tauy = tab_cntrl(29) ELSE fxyhypb = . FALSE . ysinus = . FALSE . IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE. ENDIF day_ini = tab_cntrl(30) itau_dyn = tab_cntrl(31) c ................................................................. c c c PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa cAl1 Print*,'day_ref,annee_ref,day_ini,itau_dyn', & day_ref,annee_ref,day_ini,itau_dyn c Lecture des champs c plev(1,klev+1)=0. CALL get_field("plev",plev,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("play",play,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("phi",phi,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("phis",phis,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("presnivs",presnivs,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("ucov",ucov,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("vcov",vcov,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("temp",temp,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' CALL get_field("omega2",omega2,found) IF (.NOT. found) PRINT*, modname//'Le champ est absent' Do iq=1,nqtot CALL get_field("q"//nmq(iq),q(:,:,iq),found) IF (.NOT. found) & PRINT*, modname//'Le champ est absent' EndDo CALL close_startphy print*,' close startphy' . ,fichnom,play(1,1),play(1,klev),temp(1,klev) c RETURN END ! ! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$ ! c SUBROUTINE dyn1dredem(fichnom,plev,play,phi,phis,presnivs, & ucov,vcov,temp,q,omega2) USE dimphy USE mod_grid_phy_lmdz USE mod_phys_lmdz_para USE phys_state_var_mod USE iostart USE infotrac use control_mod IMPLICIT NONE c======================================================= c Ecriture du fichier de redemarrage sous format NetCDF c======================================================= c Declarations: c ------------- #include "dimensions.h" #include "comconst.h" #include "temps.h" !!#include "control.h" #include "logic.h" #include "netcdf.inc" c Arguments: c ---------- CHARACTER*(*) fichnom REAL time cAl1 plev tronque pour .nc mais plev(klev+1):=0 real :: plev(klon,klev),play (klon,klev),phi(klon,klev) real :: presnivs(klon,klev) real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev) real :: q(klon,klev,nqtot) real :: omega2(klon,klev),rho(klon,klev+1) c real :: ug(klev),vg(klev),fcoriolis real :: phis(klon) c Variables locales pour NetCDF: c ------------------------------ INTEGER nid, nvarid INTEGER idim_s INTEGER ierr, ierr_file INTEGER iq,l INTEGER length PARAMETER (length = 100) REAL tab_cntrl(length) ! tableau des parametres du run character*4 nmq(nqtot) character*20 modname character*80 abort_message c INTEGER nb SAVE nb DATA nb / 0 / REAL zan0,zjulian,hours INTEGER yyears0,jjour0, mmois0 character*30 unites cDbg CALL open_restartphy(fichnom) print*,'redm1 ',fichnom,klon,klev,nqtot nmq(1)="vap" nmq(2)="cond" nmq(3)="tra1" nmq(4)="tra2" modname = 'dyn1dredem' ierr = NF_OPEN(fichnom, NF_WRITE, nid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Pb. d ouverture "//fichnom CALL abort ENDIF DO l=1,length tab_cntrl(l) = 0. ENDDO tab_cntrl(1) = FLOAT(iim) tab_cntrl(2) = FLOAT(jjm) tab_cntrl(3) = FLOAT(llm) tab_cntrl(4) = FLOAT(day_ref) tab_cntrl(5) = FLOAT(annee_ref) tab_cntrl(6) = rad tab_cntrl(7) = omeg tab_cntrl(8) = g tab_cntrl(9) = cpp tab_cntrl(10) = kappa tab_cntrl(11) = daysec tab_cntrl(12) = dtvr c tab_cntrl(13) = etot0 c tab_cntrl(14) = ptot0 c tab_cntrl(15) = ztot0 c tab_cntrl(16) = stot0 c tab_cntrl(17) = ang0 c tab_cntrl(18) = pa c tab_cntrl(19) = preff c c ..... parametres pour le zoom ...... c tab_cntrl(20) = clon c tab_cntrl(21) = clat c tab_cntrl(22) = grossismx c tab_cntrl(23) = grossismy c IF ( fxyhypb ) THEN tab_cntrl(24) = 1. c tab_cntrl(25) = dzoomx c tab_cntrl(26) = dzoomy tab_cntrl(27) = 0. c tab_cntrl(28) = taux c tab_cntrl(29) = tauy ELSE tab_cntrl(24) = 0. c tab_cntrl(25) = dzoomx c tab_cntrl(26) = dzoomy tab_cntrl(27) = 0. tab_cntrl(28) = 0. tab_cntrl(29) = 0. IF( ysinus ) tab_cntrl(27) = 1. ENDIF CAl1 iday_end -> day_end tab_cntrl(30) = FLOAT(day_end) tab_cntrl(31) = FLOAT(itau_dyn + itaufin) c CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl) c c Ecriture/extension de la coordonnee temps nb = nb + 1 c Ecriture des champs c CALL put_field("plev","p interfaces sauf la nulle",plev) CALL put_field("play","",play) CALL put_field("phi","geopotentielle",phi) CALL put_field("phis","geopotentiell de surface",phis) CALL put_field("presnivs","",presnivs) CALL put_field("ucov","",ucov) CALL put_field("vcov","",vcov) CALL put_field("temp","",temp) CALL put_field("omega2","",omega2) Do iq=1,nqtot CALL put_field("q"//nmq(iq),"eau vap ou condens et traceurs", . q(:,:,iq)) EndDo CALL close_restartphy c RETURN END SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn) IMPLICIT NONE !======================================================================= ! passage d'un champ de la grille scalaire a la grille physique !======================================================================= !----------------------------------------------------------------------- ! declarations: ! ------------- INTEGER im,jm,ngrid,nfield REAL pdyn(im,jm,nfield) REAL pfi(ngrid,nfield) INTEGER i,j,ifield,ig !----------------------------------------------------------------------- ! calcul: ! ------- DO ifield=1,nfield ! traitement des poles DO i=1,im pdyn(i,1,ifield)=pfi(1,ifield) pdyn(i,jm,ifield)=pfi(ngrid,ifield) ENDDO ! traitement des point normaux DO j=2,jm-1 ig=2+(j-2)*(im-1) CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1) pdyn(im,j,ifield)=pdyn(1,j,ifield) ENDDO ENDDO RETURN END SUBROUTINE abort_gcm(modname, message, ierr) USE IOIPSL ! ! Stops the simulation cleanly, closing files and printing various ! comments ! ! Input: modname = name of calling program ! message = stuff to print ! ierr = severity of situation ( = 0 normal ) character*20 modname integer ierr character*80 message write(*,*) 'in abort_gcm' call histclo ! call histclo(2) ! call histclo(3) ! call histclo(4) ! call histclo(5) write(*,*) 'out of histclo' write(*,*) 'Stopping in ', modname write(*,*) 'Reason = ',message call getin_dump ! if (ierr .eq. 0) then write(*,*) 'Everything is cool' else write(*,*) 'Houston, we have a problem ', ierr endif STOP END REAL FUNCTION fq_sat(kelvin, millibar) ! IMPLICIT none !====================================================================== ! Autheur(s): Z.X. Li (LMD/CNRS) ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.) !====================================================================== ! Arguments: ! kelvin---input-R: temperature en Kelvin ! millibar--input-R: pression en mb ! ! fq_sat----output-R: vapeur d'eau saturante en kg/kg !====================================================================== ! REAL kelvin, millibar ! REAL r2es PARAMETER (r2es=611.14 *18.0153/28.9644) ! REAL r3les, r3ies, r3es PARAMETER (R3LES=17.269) PARAMETER (R3IES=21.875) ! REAL r4les, r4ies, r4es PARAMETER (R4LES=35.86) PARAMETER (R4IES=7.66) ! REAL rtt PARAMETER (rtt=273.16) ! REAL retv PARAMETER (retv=28.9644/18.0153 - 1.0) ! REAL zqsat REAL temp, pres ! ------------------------------------------------------------------ ! ! temp = kelvin pres = millibar * 100.0 ! write(*,*)'kelvin,millibar=',kelvin,millibar ! write(*,*)'temp,pres=',temp,pres ! IF (temp .LE. rtt) THEN r3es = r3ies r4es = r4ies ELSE r3es = r3les r4es = r4les ENDIF ! zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) ) zqsat=MIN(0.5,ZQSAT) zqsat=zqsat/(1.-retv *zqsat) ! fq_sat = zqsat ! RETURN END subroutine wrgradsfi(if,nl,field,name,titlevar) implicit none ! Declarations #include "gradsdef.h" ! arguments integer if,nl real field(imx*jmx*lmx) character*10 name,file character*10 titlevar ! local integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf logical writectl writectl=.false. ! print*,if,iid(if),jid(if),ifd(if),jfd(if) iii=iid(if) iji=jid(if) iif=ifd(if) ijf=jfd(if) im=iif-iii+1 jm=ijf-iji+1 lm=lmd(if) ! print*,'im,jm,lm,name,firsttime(if)' ! print*,im,jm,lm,name,firsttime(if) if(firsttime(if)) then if(name.eq.var(1,if)) then firsttime(if)=.false. ivar(if)=1 print*,'fin de l initialiation de l ecriture du fichier' print*,file print*,'fichier no: ',if print*,'unit ',unit(if) print*,'nvar ',nvar(if) print*,'vars ',(var(iv,if),iv=1,nvar(if)) else ivar(if)=ivar(if)+1 nvar(if)=ivar(if) var(ivar(if),if)=name tvar(ivar(if),if)=trim(titlevar) nld(ivar(if),if)=nl print*,'initialisation ecriture de ',var(ivar(if),if) print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if) endif writectl=.true. itime(if)=1 else ivar(if)=mod(ivar(if),nvar(if))+1 if (ivar(if).eq.nvar(if)) then writectl=.true. itime(if)=itime(if)+1 endif if(var(ivar(if),if).ne.name) then print*,'Il faut stoker la meme succession de champs a chaque' print*,'pas de temps' print*,'fichier no: ',if print*,'unit ',unit(if) print*,'nvar ',nvar(if) print*,'vars ',(var(iv,if),iv=1,nvar(if)) stop endif endif ! print*,'ivar(if),nvar(if),var(ivar(if),if),writectl' ! print*,ivar(if),nvar(if),var(ivar(if),if),writectl do l=1,nl irec(if)=irec(if)+1 ! print*,'Ecrit rec=',irec(if),iii,iif,iji,ijf, ! s (l-1)*imd(if)*jmd(if)+(iji-1)*imd(if)+iii ! s ,(l-1)*imd(if)*jmd(if)+(ijf-1)*imd(if)+iif write(unit(if)+1,rec=irec(if)) s ((field((l-1)*imd(if)*jmd(if)+(j-1)*imd(if)+i) s ,i=iii,iif),j=iji,ijf) enddo if (writectl) then file=fichier(if) ! WARNING! on reecrase le fichier .ctl a chaque ecriture open(unit(if),file=trim(file)//'.ctl', & form='formatted',status='unknown') write(unit(if),'(a5,1x,a40)') & 'DSET ','^'//trim(file)//'.dat' write(unit(if),'(a12)') 'UNDEF 1.0E30' write(unit(if),'(a5,1x,a40)') 'TITLE ',title(if) call formcoord(unit(if),im,xd(iii,if),1.,.false.,'XDEF') call formcoord(unit(if),jm,yd(iji,if),1.,.true.,'YDEF') call formcoord(unit(if),lm,zd(1,if),1.,.false.,'ZDEF') write(unit(if),'(a4,i10,a30)') & 'TDEF ',itime(if),' LINEAR 07AUG1998 30MN ' write(unit(if),'(a4,2x,i5)') 'VARS',nvar(if) do iv=1,nvar(if) ! print*,'if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)' ! print*,if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if) write(unit(if),1000) var(iv,if),nld(iv,if)-1/nld(iv,if) & ,99,tvar(iv,if) enddo write(unit(if),'(a7)') 'ENDVARS' ! 1000 format(a5,3x,i4,i3,1x,a39) close(unit(if)) endif ! writectl return END subroutine inigrads(if,im s ,x,fx,xmin,xmax,jm,y,ymin,ymax,fy,lm,z,fz s ,dt,file,titlel) implicit none integer if,im,jm,lm,i,j,l real x(im),y(jm),z(lm),fx,fy,fz,dt real xmin,xmax,ymin,ymax integer nf character file*10,titlel*40 #include "gradsdef.h" data unit/24,32,34,36,38,40,42,44,46,48/ data nf/0/ if (if.le.nf) stop'verifier les appels a inigrads' print*,'Entree dans inigrads' nf=if title(if)=titlel ivar(if)=0 fichier(if)=trim(file) firsttime(if)=.true. dtime(if)=dt iid(if)=1 ifd(if)=im imd(if)=im do i=1,im xd(i,if)=x(i)*fx if(xd(i,if).lt.xmin) iid(if)=i+1 if(xd(i,if).le.xmax) ifd(if)=i enddo print*,'On stoke du point ',iid(if),' a ',ifd(if),' en x' jid(if)=1 jfd(if)=jm jmd(if)=jm do j=1,jm yd(j,if)=y(j)*fy if(yd(j,if).gt.ymax) jid(if)=j+1 if(yd(j,if).ge.ymin) jfd(if)=j enddo print*,'On stoke du point ',jid(if),' a ',jfd(if),' en y' print*,'Open de dat' print*,'file=',file print*,'fichier(if)=',fichier(if) print*,4*(ifd(if)-iid(if))*(jfd(if)-jid(if)) print*,trim(file)//'.dat' OPEN (unit(if)+1,FILE=trim(file)//'.dat', s FORM='UNFORMATTED', s ACCESS='DIRECT' s ,RECL=4*(ifd(if)-iid(if)+1)*(jfd(if)-jid(if)+1)) print*,'Open de dat ok' lmd(if)=lm do l=1,lm zd(l,if)=z(l)*fz enddo irec(if)=0 !CR ! print*,if,imd(if),jmd(if),lmd(if) ! print*,'if,imd(if),jmd(if),lmd(if)' return end SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi) IMPLICIT NONE !======================================================================= ! passage d'un champ de la grille scalaire a la grille physique !======================================================================= !----------------------------------------------------------------------- ! declarations: ! ------------- INTEGER im,jm,ngrid,nfield REAL pdyn(im,jm,nfield) REAL pfi(ngrid,nfield) INTEGER j,ifield,ig !----------------------------------------------------------------------- ! calcul: ! ------- IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1) s STOP 'probleme de dim' ! traitement des poles CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid) CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid) ! traitement des point normaux DO ifield=1,nfield DO j=2,jm-1 ig=2+(j-2)*(im-1) CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1) ENDDO ENDDO RETURN END SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) ! Ancienne version disvert dont on a modifie nom pour utiliser ! le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes) ! (MPL 18092012) ! ! Auteur : P. Le Van . ! IMPLICIT NONE #include "dimensions.h" #include "paramet.h" ! !======================================================================= ! ! ! s = sigma ** kappa : coordonnee verticale ! dsig(l) : epaisseur de la couche l ds la coord. s ! sig(l) : sigma a l'interface des couches l et l-1 ! ds(l) : distance entre les couches l et l-1 en coord.s ! !======================================================================= ! REAL pa,preff REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1) REAL presnivs(llm) ! ! declarations: ! ------------- ! REAL sig(llm+1),dsig(llm) ! INTEGER l REAL snorm REAL alpha,beta,gama,delta,deltaz,h INTEGER np,ierr REAL pi,x !----------------------------------------------------------------------- ! pi=2.*ASIN(1.) OPEN(99,file='sigma.def',status='old',form='formatted', s iostat=ierr) !----------------------------------------------------------------------- ! cas 1 on lit les options dans sigma.def: ! ---------------------------------------- IF (ierr.eq.0) THEN print*,'WARNING!!! on lit les options dans sigma.def' READ(99,*) deltaz READ(99,*) h READ(99,*) beta READ(99,*) gama READ(99,*) delta READ(99,*) np CLOSE(99) alpha=deltaz/(llm*h) ! DO 1 l = 1, llm dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))* $ ( (tanh(gama*l)/tanh(gama*llm))**np + $ (1.-l/FLOAT(llm))*delta ) 1 CONTINUE sig(1)=1. DO 101 l=1,llm-1 sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l)) 101 CONTINUE sig(llm+1)=0. DO 2 l = 1, llm dsig(l) = sig(l)-sig(l+1) 2 CONTINUE ! ELSE !----------------------------------------------------------------------- ! cas 2 ancienne discretisation (LMD5...): ! ---------------------------------------- PRINT*,'WARNING!!! Ancienne discretisation verticale' h=7. snorm = 0. DO l = 1, llm x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1) dsig(l) = 1.0 + 7.0 * SIN(x)**2 snorm = snorm + dsig(l) ENDDO snorm = 1./snorm DO l = 1, llm dsig(l) = dsig(l)*snorm ENDDO sig(llm+1) = 0. DO l = llm, 1, -1 sig(l) = sig(l+1) + dsig(l) ENDDO ENDIF DO l=1,llm nivsigs(l) = FLOAT(l) ENDDO DO l=1,llmp1 nivsig(l)= FLOAT(l) ENDDO ! ! .... Calculs de ap(l) et de bp(l) .... ! ......................................... ! ! ! ..... pa et preff sont lus sur les fichiers start par lectba ..... ! bp(llmp1) = 0. DO l = 1, llm !c !cc ap(l) = 0. !cc bp(l) = sig(l) bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) ap(l) = pa * ( sig(l) - bp(l) ) ! ENDDO ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) PRINT *,' BP ' PRINT *, bp PRINT *,' AP ' PRINT *, ap DO l = 1, llm dpres(l) = bp(l) - bp(l+1) presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) ENDDO PRINT *,' PRESNIVS ' PRINT *,presnivs RETURN END !====================================================================== SUBROUTINE read_tsurf1d(knon,knindex,sst_out) ! This subroutine specifies the surface temperature to be used in 1D simulations USE dimphy, ONLY : klon INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model INTEGER :: i ! COMMON defined in lmdz1d.F: real ts_cur common /sst_forcing/ts_cur DO i = 1, knon sst_out(i) = ts_cur ENDDO END SUBROUTINE read_tsurf1d !=============================================================== subroutine advect_vert(llm,w,dt,q,plev) !=============================================================== ! Schema amont pour l'advection verticale en 1D ! w est la vitesse verticale dp/dt en Pa/s ! Traitement en volumes finis ! d / dt ( zm q ) = delta_z ( omega q ) ! d / dt ( zm ) = delta_z ( omega ) ! avec zm = delta_z ( p ) ! si * designe la valeur au pas de temps t+dt ! zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l) ! zm*(l) -zm(l) = w(l+1) - w(l) ! avec w=omega * dt !--------------------------------------------------------------- implicit none ! arguments integer llm real w(llm+1),q(llm),plev(llm+1),dt ! local integer l real zwq(llm+1),zm(llm+1),zw(llm+1) real qold !--------------------------------------------------------------- do l=1,llm zw(l)=dt*w(l) zm(l)=plev(l)-plev(l+1) zwq(l)=q(l)*zw(l) enddo zwq(llm+1)=0. zw(llm+1)=0. do l=1,llm qold=q(l) q(l)=(q(l)*zm(l)+zwq(l+1)-zwq(l))/(zm(l)+zw(l+1)-zw(l)) print*,'ADV Q ',zm(l),zw(l),zwq(l),qold,q(l) enddo return end !=============================================================== SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va, ! q,temp,u,v, ! play,plev) !itlmd !---------------------------------------------------------------------- ! Calcul de l'advection verticale (ascendance et subsidence) de ! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur ! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou ! sans WTG rajouter une advection horizontale !---------------------------------------------------------------------- implicit none #include "YOMCST.h" ! argument integer llm real omega(llm+1),d_t_va(llm), d_q_va(llm,3) real d_u_va(llm), d_v_va(llm) real q(llm,3),temp(llm) real u(llm),v(llm) real play(llm),plev(llm+1) ! interne integer l real alpha,omgdown,omgup do l= 1,llm if(l.eq.1) then !si omgup pour la couche 1, alors tendance nulle omgdown=max(omega(2),0.0) alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)* & (1.+q(l,1))) d_t_va(l)= alpha*(omgdown)- & omgdown*(temp(l)-temp(l+1)) & /(play(l)-play(l+1)) d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:)) & /(play(l)-play(l+1)) d_u_va(l)= -omgdown*(u(l)-u(l+1)) & /(play(l)-play(l+1)) d_v_va(l)= -omgdown*(v(l)-v(l+1)) & /(play(l)-play(l+1)) elseif(l.eq.llm) then omgup=min(omega(l),0.0) alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)* & (1.+q(l,1))) d_t_va(l)= alpha*(omgup)- !bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l)) & omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l)) d_q_va(l,:)= -omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l)) d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l)) d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l)) else omgup=min(omega(l),0.0) omgdown=max(omega(l+1),0.0) alpha = rkappa*temp(l)*(1.+q(l,1)*rv/rd)/(play(l)* & (1.+q(l,1))) d_t_va(l)= alpha*(omgup+omgdown)- & omgdown*(temp(l)-temp(l+1)) & /(play(l)-play(l+1))- !bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l)) & omgup*(temp(l-1)-temp(l))/(play(l-1)-play(l)) ! print*, ' ??? ' d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:)) & /(play(l)-play(l+1))- & omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l)) d_u_va(l)= -omgdown*(u(l)-u(l+1)) & /(play(l)-play(l+1))- & omgup*(u(l-1)-u(l))/(play(l-1)-play(l)) d_v_va(l)= -omgdown*(v(l)-v(l+1)) & /(play(l)-play(l+1))- & omgup*(v(l-1)-v(l))/(play(l-1)-play(l)) endif enddo !fin itlmd return end ! SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va, SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va, ! q,temp,u,v,play) !itlmd !---------------------------------------------------------------------- ! Calcul de l'advection verticale (ascendance et subsidence) de ! température et d'humidité. Hypothèse : ce qui rentre de l'extérieur ! a les mêmes caractéristiques que l'air de la colonne 1D (WTG) ou ! sans WTG rajouter une advection horizontale !---------------------------------------------------------------------- implicit none #include "YOMCST.h" ! argument integer llm,nqtot real omega(llm+1),d_t_va(llm), d_q_va(llm,nqtot) ! real d_u_va(llm), d_v_va(llm) real q(llm,nqtot),temp(llm) real u(llm),v(llm) real play(llm) real cor(llm) ! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm) real dph(llm),dqdp(llm),dtdp(llm) ! interne integer l,k real alpha,omdn,omup ! dudp=0. ! dvdp=0. dqdp=0. dtdp=0. ! d_u_va=0. ! d_v_va=0. cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1))) do k=2,llm-1 dph (k-1) = (play(k )- play(k-1 )) ! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1) ! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1) dqdp (k-1) = (q (k,1)- q (k-1,1))/dph(k-1) dtdp (k-1) = (temp(k )- temp(k-1 ))/dph(k-1) enddo ! dudp ( llm ) = dudp ( llm-1 ) ! dvdp ( llm ) = dvdp ( llm-1 ) dqdp ( llm ) = dqdp ( llm-1 ) dtdp ( llm ) = dtdp ( llm-1 ) do k=2,llm-1 omdn=max(0.0,omega(k+1)) omup=min(0.0,omega( k )) ! d_u_va(k) = -omdn*dudp(k)-omup*dudp(k-1) ! d_v_va(k) = -omdn*dvdp(k)-omup*dvdp(k-1) d_q_va(k,1)= -omdn*dqdp(k)-omup*dqdp(k-1) d_t_va(k) = -omdn*dtdp(k)-omup*dtdp(k-1) : +(omup+omdn)*cor(k) enddo omdn=max(0.0,omega( 2 )) omup=min(0.0,omega(llm)) ! d_u_va( 1 ) = -omdn*dudp( 1 ) ! d_u_va(llm) = -omup*dudp(llm) ! d_v_va( 1 ) = -omdn*dvdp( 1 ) ! d_v_va(llm) = -omup*dvdp(llm) d_q_va( 1 ,1) = -omdn*dqdp( 1 ) d_q_va(llm,1) = -omup*dqdp(llm) d_t_va( 1 ) = -omdn*dtdp( 1 )+omdn*cor( 1 ) d_t_va(llm) = -omup*dtdp(llm)!+omup*cor(llm) ! if(abs(rlat(1))>10.) then ! Calculate the tendency due agestrophic motions ! du_age = fcoriolis*(v-vg) ! dv_age = fcoriolis*(ug-u) ! endif ! call writefield_phy('d_t_va',d_t_va,llm) return end !====================================================================== SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga : ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga : ,ht_toga,vt_toga,hq_toga,vq_toga) implicit none c------------------------------------------------------------------------- c Read TOGA-COARE forcing data c------------------------------------------------------------------------- integer nlev_toga,nt_toga real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga) real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga) real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga) real w_toga(nlev_toga,nt_toga) real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga) real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga) character*80 fich_toga integer no,l,k,ip real riy,rim,rid,rih,bid integer iy,im,id,ih real plev_min plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa open(21,file=trim(fich_toga),form='formatted') read(21,'(a)') do ip = 1, nt_toga read(21,'(a)') read(21,'(a)') read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid read(21,'(a)') read(21,'(a)') do k = 1, nlev_toga read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip) : ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip) : ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip) ! conversion in SI units: t_toga(k,ip)=t_toga(k,ip)+273.15 ! K q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s ! no water vapour tendency above 55 hPa if (plev_toga(k,ip) .lt. plev_min) then q_toga(k,ip) = 0. hq_toga(k,ip) = 0. vq_toga(k,ip) =0. endif enddo ts_toga(ip)=ts_toga(ip)+273.15 ! K enddo close(21) 223 format(4i3,6f8.2) 226 format(f7.1,1x,10f8.2) 227 format(f7.1,1x,1p,4e11.3) 230 format(6f9.3,4e11.3) return end c------------------------------------------------------------------------- SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu) implicit none c------------------------------------------------------------------------- c Read I.SANDU case forcing data c------------------------------------------------------------------------- integer nlev_sandu,nt_sandu real ts_sandu(nt_sandu) character*80 fich_sandu integer no,l,k,ip real riy,rim,rid,rih,bid integer iy,im,id,ih real plev_min plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa open(21,file=trim(fich_sandu),form='formatted') read(21,'(a)') do ip = 1, nt_sandu read(21,'(a)') read(21,'(a)') read(21,223) iy, im, id, ih, ts_sandu(ip) print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip) enddo close(21) 223 format(4i3,f8.2) 226 format(f7.1,1x,10f8.2) 227 format(f7.1,1x,1p,4e11.3) 230 format(6f9.3,4e11.3) return end !===================================================================== c------------------------------------------------------------------------- SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex, : ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex) implicit none c------------------------------------------------------------------------- c Read Astex case forcing data c------------------------------------------------------------------------- integer nlev_astex,nt_astex real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex) real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex) character*80 fich_astex integer no,l,k,ip real riy,rim,rid,rih,bid integer iy,im,id,ih real plev_min plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa open(21,file=trim(fich_astex),form='formatted') read(21,'(a)') read(21,'(a)') do ip = 1, nt_astex read(21,'(a)') read(21,'(a)') read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip), :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip) ts_astex(ip)=ts_astex(ip)+273.15 print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip), :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip) enddo close(21) 223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2) 226 format(f7.1,1x,10f8.2) 227 format(f7.1,1x,1p,4e11.3) 230 format(6f9.3,4e11.3) return end !===================================================================== subroutine read_twpice(fich_twpice,nlevel,ntime : ,T_srf,plev,T,q,u,v,omega : ,T_adv_h,T_adv_v,q_adv_h,q_adv_v) !program reading forcings of the TWP-ICE experiment ! use netcdf implicit none #include "netcdf.inc" integer ntime,nlevel integer l,k character*80 :: fich_twpice real*8 time(ntime) real*8 lat, lon, alt, phis real*8 lev(nlevel) real*8 plev(nlevel,ntime) real*8 T(nlevel,ntime) real*8 q(nlevel,ntime),u(nlevel,ntime) real*8 v(nlevel,ntime) real*8 omega(nlevel,ntime), div(nlevel,ntime) real*8 T_adv_h(nlevel,ntime) real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime) real*8 q_adv_v(nlevel,ntime) real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime) real*8 s_adv_v(nlevel,ntime) real*8 p_srf_aver(ntime), p_srf_center(ntime) real*8 T_srf(ntime) integer nid, ierr integer nbvar3d parameter(nbvar3d=20) integer var3didin(nbvar3d) ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid) if (ierr.NE.NF_NOERR) then write(*,*) 'ERROR: Pb opening forcings cdf file ' write(*,*) NF_STRERROR(ierr) stop "" endif ierr=NF_INQ_VARID(nid,"lat",var3didin(1)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'lat' endif ierr=NF_INQ_VARID(nid,"lon",var3didin(2)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'lon' endif ierr=NF_INQ_VARID(nid,"alt",var3didin(3)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'alt' endif ierr=NF_INQ_VARID(nid,"phis",var3didin(4)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'phis' endif ierr=NF_INQ_VARID(nid,"T",var3didin(5)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'T' endif ierr=NF_INQ_VARID(nid,"q",var3didin(6)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'q' endif ierr=NF_INQ_VARID(nid,"u",var3didin(7)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'u' endif ierr=NF_INQ_VARID(nid,"v",var3didin(8)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'v' endif ierr=NF_INQ_VARID(nid,"omega",var3didin(9)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'omega' endif ierr=NF_INQ_VARID(nid,"div",var3didin(10)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'div' endif ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'T_adv_h' endif ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'T_adv_v' endif ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'q_adv_h' endif ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'q_adv_v' endif ierr=NF_INQ_VARID(nid,"s",var3didin(15)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 's' endif ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 's_adv_h' endif ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 's_adv_v' endif ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'p_srf_aver' endif ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'p_srf_center' endif ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20)) if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop 'T_srf' endif !dimensions lecture call catchaxis(nid,ntime,nlevel,time,lev,ierr) !pressure do l=1,ntime do k=1,nlevel plev(k,l)=lev(k) enddo enddo #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat) #else ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture lat ok',lat #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon) #else ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture lon ok',lon #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt) #else ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture alt ok',alt #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis) #else ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture phis ok',phis #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T) #else ierr = NF_GET_VAR_REAL(nid,var3didin(5),T) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture T ok' #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q) #else ierr = NF_GET_VAR_REAL(nid,var3didin(6),q) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture q ok' !q in kg/kg do l=1,ntime do k=1,nlevel q(k,l)=q(k,l)/1000. enddo enddo #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u) #else ierr = NF_GET_VAR_REAL(nid,var3didin(7),u) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture u ok' #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v) #else ierr = NF_GET_VAR_REAL(nid,var3didin(8),v) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture v ok' #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega) #else ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture omega ok' !omega in mb/hour do l=1,ntime do k=1,nlevel omega(k,l)=omega(k,l)*100./3600. enddo enddo #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div) #else ierr = NF_GET_VAR_REAL(nid,var3didin(10),div) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture div ok' #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h) #else ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture T_adv_h ok' !T adv in K/s do l=1,ntime do k=1,nlevel T_adv_h(k,l)=T_adv_h(k,l)/3600. enddo enddo #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v) #else ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture T_adv_v ok' !T adv in K/s do l=1,ntime do k=1,nlevel T_adv_v(k,l)=T_adv_v(k,l)/3600. enddo enddo #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h) #else ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture q_adv_h ok' !q adv in kg/kg/s do l=1,ntime do k=1,nlevel q_adv_h(k,l)=q_adv_h(k,l)/1000./3600. enddo enddo #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v) #else ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture q_adv_v ok' !q adv in kg/kg/s do l=1,ntime do k=1,nlevel q_adv_v(k,l)=q_adv_v(k,l)/1000./3600. enddo enddo #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s) #else ierr = NF_GET_VAR_REAL(nid,var3didin(15),s) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h) #else ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v) #else ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver) #else ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center) #else ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf) #else ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf) #endif if(ierr/=NF_NOERR) then write(*,*) NF_STRERROR(ierr) stop "getvarup" endif ! write(*,*)'lecture T_srf ok', T_srf return end subroutine read_twpice !===================================================================== subroutine catchaxis(nid,ttm,llm,time,lev,ierr) ! use netcdf implicit none #include "netcdf.inc" integer nid,ttm,llm real*8 time(ttm) real*8 lev(llm) integer ierr integer i integer timevar,levvar integer timelen,levlen integer timedimin,levdimin ! Control & lecture on dimensions ! =============================== ierr=NF_INQ_DIMID(nid,"time",timedimin) ierr=NF_INQ_VARID(nid,"time",timevar) if (ierr.NE.NF_NOERR) then write(*,*) 'ERROR: Field