MODULE lmdz_1dutils IMPLICIT NONE; PRIVATE PUBLIC fq_sat, conf_unicol, dyn1deta0, dyn1dredem, & disvert0, advect_vert, advect_va, lstendh, nudge_rht_init, nudge_uv_init, & nudge_rht, nudge_uv, interp2_case_vertical CONTAINS 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, INTENT(IN) :: 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 <= 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 END FUNCTION fq_sat SUBROUTINE conf_unicol USE IOIPSL USE lmdz_print_control, ONLY: lunout USE lmdz_flux_arp, ONLY: fsens, flat, betaevap, ust, tg, ok_flux_surf, ok_prescr_ust, ok_prescr_beta, ok_forc_tsurf USE lmdz_fcs_gcssold, ONLY: imp_fcg_gcssold, ts_fcg_gcssold, Tp_fcg_gcssold, Tp_ini_gcssold, xTurb_fcg_gcssold USE lmdz_tsoilnudge, ONLY: nudge_tsoil, isoil_nudge, Tsoil_nudge, tau_soil_nudge USE lmdz_compar1d IMPLICIT NONE !----------------------------------------------------------------------- ! Auteurs : A. Lahellec . ! ------------------------------------------------------------------- ! ......... Initilisation parametres du lmdz1D .......... !--------------------------------------------------------------------- ! initialisations: ! ---------------- !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 debogage !Config Def = 0 !Config Help = Niveau d'impression pour le debogage !Config (0 = minimum d'impression) ! prt_level = 0 ! CALL getin('prt_level',prt_level) !----------------------------------------------------------------------- ! Parametres de controle du run: !----------------------------------------------------------------------- !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. ! = 10 ==> forcing_case = .TRUE. ! initial profiles from case.nc file ! = 40 ==> forcing_GCSSold = .TRUE. ! initial profile from GCSS file ! LS convergence imposed from GCSS file ! = 50 ==> forcing_fire = .TRUE. ! = 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 ! > 100 ==> forcing_case = .TRUE. or forcing_case2 = .TRUE. ! initial profiles from case.nc file 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 ==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 !Parametres de forcage !Config Key = tend_t !Config Desc = forcage ou non par advection de T !Config Def = false !Config Help = forcage ou non par advection de T tend_t = 0 CALL getin('tend_t', tend_t) !Config Key = tend_q !Config Desc = forcage ou non par advection de q !Config Def = false !Config Help = forcage ou non par advection de q tend_q = 0 CALL getin('tend_q', tend_q) !Config Key = tend_u !Config Desc = forcage ou non par advection de u !Config Def = false !Config Help = forcage ou non par advection de u tend_u = 0 CALL getin('tend_u', tend_u) !Config Key = tend_v !Config Desc = forcage ou non par advection de v !Config Def = false !Config Help = forcage ou non par advection de v tend_v = 0 CALL getin('tend_v', tend_v) !Config Key = tend_w !Config Desc = forcage ou non par vitesse verticale !Config Def = false !Config Help = forcage ou non par vitesse verticale tend_w = 0 CALL getin('tend_w', tend_w) !Config Key = tend_rayo !Config Desc = forcage ou non par dtrad !Config Def = false !Config Help = forcage ou non par dtrad tend_rayo = 0 CALL getin('tend_rayo', tend_rayo) !Config Key = nudge_t !Config Desc = constante de nudging de T !Config Def = false !Config Help = constante de nudging de T nudge_t = 0. CALL getin('nudge_t', nudge_t) !Config Key = nudge_q !Config Desc = constante de nudging de q !Config Def = false !Config Help = constante de nudging de q nudge_q = 0. CALL getin('nudge_q', nudge_q) !Config Key = nudge_u !Config Desc = constante de nudging de u !Config Def = false !Config Help = constante de nudging de u nudge_u = 0. CALL getin('nudge_u', nudge_u) !Config Key = nudge_v !Config Desc = constante de nudging de v !Config Def = false !Config Help = constante de nudging de v nudge_v = 0. CALL getin('nudge_v', nudge_v) !Config Key = nudge_w !Config Desc = constante de nudging de w !Config Def = false !Config Help = constante de nudging de w nudge_w = 0. CALL getin('nudge_w', nudge_w) !Config Key = iflag_nudge !Config Desc = atmospheric nudging ttype (decimal code) !Config Def = 0 !Config Help = 0 ==> no nudging ! If digit number n of iflag_nudge is set, then nudging of type n is on ! If digit number n of iflag_nudge is not set, then nudging of type n is off ! (digits are numbered from the right) iflag_nudge = 0 CALL getin('iflag_nudge', iflag_nudge) !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_forc_tsurf !Config Desc = forcage ou non par la Ts !Config Def = false !Config Help = forcage ou non par la Ts ok_forc_tsurf = .FALSE. CALL getin('ok_forc_tsurf', ok_forc_tsurf) !Config Key = ok_prescr_ust !Config Desc = ustar impose ou non !Config Def = false !Config Help = ustar impose ou non ok_prescr_ust = .FALSE. CALL getin('ok_prescr_ust', ok_prescr_ust) !Config Key = ok_prescr_beta !Config Desc = betaevap impose ou non !Config Def = false !Config Help = betaevap impose ou non ok_prescr_beta = .FALSE. CALL getin('ok_prescr_beta', ok_prescr_beta) !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 = time_ini = 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 = surface temperature 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) ! EV pour accord avec format standard CALL getin('zorog', 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) ! FH/2020/04/08/confinement: Pour le nouveau format standard, la rugosite s'appelle z0 CALL getin('z0', rugos) !Config Key = rugosh !Config Desc = coefficient de frottement !Config Def = rugos !Config Help = calcul du Cdrag rugosh = rugos CALL getin('rugosh', rugosh) !Config Key = snowmass !Config Desc = mass de neige de la surface en kg/m2 !Config Def = 0.0000 !Config Help = snowmass snowmass = 0.0000 CALL getin('snowmass', snowmass) !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 = betaevap !Config Desc = beta for actual evaporation when prescribed !Config Def = 1.0 !Config Help = betaevap = 1. CALL getin('betaevap', betaevap) !Config Key = zpicinp !Config Desc = denivellation orographie !Config Def = 0. !Config Help = input brise zpicinp = 0. 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) !---------------------------------------------------------- ! Param??tres de for??age pour les forcages communs: ! Pour les forcages communs: ces entiers valent 0 ou 1 ! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale ! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale ! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv) ! forcages en omega, w, vent geostrophique ou ustar ! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging !---------------------------------------------------------- !Config Key = tadv !Config Desc = forcage ou non par advection totale de T !Config Def = false !Config Help = forcage ou non par advection totale de T tadv = 0 CALL getin('tadv', tadv) !Config Key = tadvv !Config Desc = forcage ou non par advection verticale de T !Config Def = false !Config Help = forcage ou non par advection verticale de T tadvv = 0 CALL getin('tadvv', tadvv) !Config Key = tadvh !Config Desc = forcage ou non par advection horizontale de T !Config Def = false !Config Help = forcage ou non par advection horizontale de T tadvh = 0 CALL getin('tadvh', tadvh) !Config Key = thadv !Config Desc = forcage ou non par advection totale de Theta !Config Def = false !Config Help = forcage ou non par advection totale de Theta thadv = 0 CALL getin('thadv', thadv) !Config Key = thadvv !Config Desc = forcage ou non par advection verticale de Theta !Config Def = false !Config Help = forcage ou non par advection verticale de Theta thadvv = 0 CALL getin('thadvv', thadvv) !Config Key = thadvh !Config Desc = forcage ou non par advection horizontale de Theta !Config Def = false !Config Help = forcage ou non par advection horizontale de Theta thadvh = 0 CALL getin('thadvh', thadvh) !Config Key = qadv !Config Desc = forcage ou non par advection totale de Q !Config Def = false !Config Help = forcage ou non par advection totale de Q qadv = 0 CALL getin('qadv', qadv) !Config Key = qadvv !Config Desc = forcage ou non par advection verticale de Q !Config Def = false !Config Help = forcage ou non par advection verticale de Q qadvv = 0 CALL getin('qadvv', qadvv) !Config Key = qadvh !Config Desc = forcage ou non par advection horizontale de Q !Config Def = false !Config Help = forcage ou non par advection horizontale de Q qadvh = 0 CALL getin('qadvh', qadvh) !Config Key = trad !Config Desc = forcage ou non par tendance radiative !Config Def = false !Config Help = forcage ou non par tendance radiative trad = 0 CALL getin('trad', trad) !Config Key = forc_omega !Config Desc = forcage ou non par omega !Config Def = false !Config Help = forcage ou non par omega forc_omega = 0 CALL getin('forc_omega', forc_omega) !Config Key = forc_u !Config Desc = forcage ou non par u !Config Def = false !Config Help = forcage ou non par u forc_u = 0 CALL getin('forc_u', forc_u) !Config Key = forc_v !Config Desc = forcage ou non par v !Config Def = false !Config Help = forcage ou non par v forc_v = 0 CALL getin('forc_v', forc_v) !Config Key = forc_w !Config Desc = forcage ou non par w !Config Def = false !Config Help = forcage ou non par w forc_w = 0 CALL getin('forc_w', forc_w) !Config Key = forc_geo !Config Desc = forcage ou non par geo !Config Def = false !Config Help = forcage ou non par geo forc_geo = 0 CALL getin('forc_geo', forc_geo) ! Meme chose que ok_precr_ust !Config Key = forc_ustar !Config Desc = forcage ou non par ustar !Config Def = false !Config Help = forcage ou non par ustar forc_ustar = 0 CALL getin('forc_ustar', forc_ustar) IF (forc_ustar == 1) ok_prescr_ust = .TRUE. !Config Key = nudging_u !Config Desc = forcage ou non par nudging sur u !Config Def = false !Config Help = forcage ou non par nudging sur u nudging_u = 0 CALL getin('nudging_u', nudging_u) !Config Key = nudging_v !Config Desc = forcage ou non par nudging sur v !Config Def = false !Config Help = forcage ou non par nudging sur v nudging_v = 0 CALL getin('nudging_v', nudging_v) !Config Key = nudging_w !Config Desc = forcage ou non par nudging sur w !Config Def = false !Config Help = forcage ou non par nudging sur w nudging_w = 0 CALL getin('nudging_w', nudging_w) ! RELIQUE ANCIENS FORMAT. ECRASE PAR LE SUIVANT !Config Key = nudging_q !Config Desc = forcage ou non par nudging sur q !Config Def = false !Config Help = forcage ou non par nudging sur q nudging_qv = 0 CALL getin('nudging_q', nudging_qv) CALL getin('nudging_qv', nudging_qv) p_nudging_u = 11000. p_nudging_v = 11000. p_nudging_t = 11000. p_nudging_qv = 11000. CALL getin('p_nudging_u', p_nudging_u) CALL getin('p_nudging_v', p_nudging_v) CALL getin('p_nudging_t', p_nudging_t) CALL getin('p_nudging_qv', p_nudging_qv) !Config Key = nudging_t !Config Desc = forcage ou non par nudging sur t !Config Def = false !Config Help = forcage ou non par nudging sur t nudging_t = 0 CALL getin('nudging_t', nudging_t) 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, *)' snowmass=', snowmass 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 WRITE(lunout, *)' tadv = ', tadv WRITE(lunout, *)' tadvv = ', tadvv WRITE(lunout, *)' tadvh = ', tadvh WRITE(lunout, *)' thadv = ', thadv WRITE(lunout, *)' thadvv = ', thadvv WRITE(lunout, *)' thadvh = ', thadvh WRITE(lunout, *)' qadv = ', qadv WRITE(lunout, *)' qadvv = ', qadvv WRITE(lunout, *)' qadvh = ', qadvh WRITE(lunout, *)' trad = ', trad WRITE(lunout, *)' forc_omega = ', forc_omega WRITE(lunout, *)' forc_w = ', forc_w WRITE(lunout, *)' forc_geo = ', forc_geo WRITE(lunout, *)' forc_ustar = ', forc_ustar WRITE(lunout, *)' nudging_u = ', nudging_u WRITE(lunout, *)' nudging_v = ', nudging_v WRITE(lunout, *)' nudging_t = ', nudging_t WRITE(lunout, *)' nudging_qv = ', nudging_qv IF (forcing_type ==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, *) END SUBROUTINE conf_unicol SUBROUTINE dyn1deta0(fichnom, plev, play, phi, phis, presnivs, & & ucov, vcov, temp, q, omega2) USE dimphy USE lmdz_grid_phy USE lmdz_phys_para USE iophy USE phys_state_var_mod USE iostart USE lmdz_writefield_phy USE lmdz_infotrac USE control_mod USE comconst_mod, ONLY: im, jm, lllm USE logic_mod, ONLY: fxyhypb, ysinus USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm IMPLICIT NONE !======================================================= ! Ecriture du fichier de redemarrage sous format NetCDF !======================================================= ! Declarations: ! ------------- ! Arguments: ! ---------- CHARACTER*(*) fichnom !Al1 plev tronque pour .nc mais plev(klev+1):=0 REAL :: plev(klon, klev + 1), 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) ! REAL :: ug(klev),vg(klev),fcoriolis REAL :: phis(klon) ! Variables locales pour NetCDF: ! ------------------------------ 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 modname = 'dyn1deta0 : ' !! nmq(1)="vap" !! nmq(2)="cond" !! do iq=3,nqtot !! WRITE(nmq(iq),'("tra",i1)') iq-2 !! enddo DO iq = 1, nqtot nmq(iq) = trim(tracers(iq)%name) ENDDO PRINT*, 'in dyn1deta0 ', fichnom, klon, klev, nqtot CALL open_startphy(fichnom) PRINT*, 'after open startphy ', fichnom, nmq ! Lecture des parametres de controle: 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) ! rad = tab_cntrl(6) ! omeg = tab_cntrl(7) ! g = tab_cntrl(8) ! cpp = tab_cntrl(9) ! kappa = tab_cntrl(10) ! daysec = tab_cntrl(11) ! dtvr = tab_cntrl(12) ! etot0 = tab_cntrl(13) ! ptot0 = tab_cntrl(14) ! ztot0 = tab_cntrl(15) ! stot0 = tab_cntrl(16) ! ang0 = tab_cntrl(17) ! pa = tab_cntrl(18) ! preff = tab_cntrl(19) ! clon = tab_cntrl(20) ! clat = tab_cntrl(21) ! grossismx = tab_cntrl(22) ! grossismy = tab_cntrl(23) IF (tab_cntrl(24)==1.) THEN fxyhypb = .TRUE. ! dzoomx = tab_cntrl(25) ! dzoomy = tab_cntrl(26) ! taux = tab_cntrl(28) ! tauy = tab_cntrl(29) ELSE fxyhypb = .FALSE. ysinus = .FALSE. IF(tab_cntrl(27)==1.) ysinus = .TRUE. ENDIF day_ini = tab_cntrl(30) itau_dyn = tab_cntrl(31) Print*, 'day_ref,annee_ref,day_ini,itau_dyn', day_ref, annee_ref, day_ini, itau_dyn ! Lecture des champs 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' plev(1, klev + 1) = 0. CALL get_field("plev", plev(:, 1:klev), 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) END SUBROUTINE dyn1deta0 SUBROUTINE dyn1dredem(fichnom, plev, play, phi, phis, presnivs, & & ucov, vcov, temp, q, omega2) USE dimphy USE lmdz_grid_phy USE lmdz_phys_para USE phys_state_var_mod USE iostart USE lmdz_infotrac USE control_mod USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad USE logic_mod, ONLY: fxyhypb, ysinus USE temps_mod, ONLY: annee_ref, day_end, day_ref, itau_dyn, itaufin USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm IMPLICIT NONE !======================================================= ! Ecriture du fichier de redemarrage sous format NetCDF !======================================================= ! Declarations: ! ------------- ! Arguments: ! ---------- CHARACTER*(*) fichnom !Al1 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) ! REAL :: ug(klev),vg(klev),fcoriolis REAL :: phis(klon) ! Variables locales pour NetCDF: ! ------------------------------ INTEGER nid INTEGER ierr 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 INTEGER pass CALL open_restartphy(fichnom) PRINT*, 'redm1 ', fichnom, klon, klev, nqtot !! nmq(1)="vap" !! nmq(2)="cond" !! nmq(3)="tra1" !! nmq(4)="tra2" DO iq = 1, nqtot nmq(iq) = trim(tracers(iq)%name) ENDDO ! modname = 'dyn1dredem' ! ierr = nf90_open(fichnom, nf90_write, nid) ! IF (ierr .NE. nf90_noerr) THEN ! abort_message="Pb. d ouverture "//fichnom ! CALL abort_gcm('Modele 1D',abort_message,1) ! 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 ! tab_cntrl(13) = etot0 ! tab_cntrl(14) = ptot0 ! tab_cntrl(15) = ztot0 ! tab_cntrl(16) = stot0 ! tab_cntrl(17) = ang0 ! tab_cntrl(18) = pa ! tab_cntrl(19) = preff ! ..... parametres pour le zoom ...... ! tab_cntrl(20) = clon ! tab_cntrl(21) = clat ! tab_cntrl(22) = grossismx ! tab_cntrl(23) = grossismy IF (fxyhypb) THEN tab_cntrl(24) = 1. ! tab_cntrl(25) = dzoomx ! tab_cntrl(26) = dzoomy tab_cntrl(27) = 0. ! tab_cntrl(28) = taux ! tab_cntrl(29) = tauy ELSE tab_cntrl(24) = 0. ! tab_cntrl(25) = dzoomx ! tab_cntrl(26) = dzoomy tab_cntrl(27) = 0. tab_cntrl(28) = 0. tab_cntrl(29) = 0. IF(ysinus) tab_cntrl(27) = 1. ENDIF !Al1 iday_end -> day_end tab_cntrl(30) = FLOAT(day_end) tab_cntrl(31) = FLOAT(itau_dyn + itaufin) DO pass = 1, 2 CALL put_var(pass, "controle", "Param. de controle Dyn1D", tab_cntrl) ! Ecriture/extension de la coordonnee temps ! Ecriture des champs CALL put_field(pass, "plev", "p interfaces sauf la nulle", plev) CALL put_field(pass, "play", "", play) CALL put_field(pass, "phi", "geopotentielle", phi) CALL put_field(pass, "phis", "geopotentiell de surface", phis) CALL put_field(pass, "presnivs", "", presnivs) CALL put_field(pass, "ucov", "", ucov) CALL put_field(pass, "vcov", "", vcov) CALL put_field(pass, "temp", "", temp) CALL put_field(pass, "omega2", "", omega2) Do iq = 1, nqtot CALL put_field(pass, "q" // nmq(iq), "eau vap ou condens et traceurs", & & q(:, :, iq)) EndDo IF (pass==1) CALL enddef_restartphy IF (pass==2) CALL close_restartphy ENDDO END SUBROUTINE dyn1dredem 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 . USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE !======================================================================= ! 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', & & iostat = ierr) !----------------------------------------------------------------------- ! cas 1 on lit les options dans sigma.def: ! ---------------------------------------- IF (ierr==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 l = 1, llm dsig(l) = (alpha + (1. - alpha) * exp(-beta * (llm - l))) * & & ((tanh(gama * l) / tanh(gama * llm))**np + & & (1. - l / FLOAT(llm)) * delta) END DO sig(1) = 1. DO l = 1, llm - 1 sig(l + 1) = sig(l) * (1. - dsig(l)) / (1. + dsig(l)) END DO sig(llm + 1) = 0. DO l = 1, llm dsig(l) = sig(l) - sig(l + 1) END DO 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 END SUBROUTINE disvert0 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 END SUBROUTINE advect_vert SUBROUTINE advect_va(llm, omega, d_t_va, d_q_va, d_u_va, d_v_va, & & q, temp, u, v, play) !itlmd !---------------------------------------------------------------------- ! Calcul de l'advection verticale (ascendance et subsidence) de ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou ! sans WTG rajouter une advection horizontale !---------------------------------------------------------------------- USE lmdz_yomcst IMPLICIT NONE ! 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) ! interne INTEGER l REAL alpha, omgdown, omgup DO l = 1, llm IF(l==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==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 END SUBROUTINE advect_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 ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou ! sans WTG rajouter une advection horizontale !---------------------------------------------------------------------- USE lmdz_yomcst IMPLICIT NONE ! 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 k REAL 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) END SUBROUTINE lstendH SubROUTINE Nudge_RHT_init(paprs, pplay, t, q, t_targ, rh_targ) ! ======================================================== USE dimphy USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE INCLUDE "FCTTRE.h" ! ======================================================== REAL paprs(klon, klevp1) REAL pplay(klon, klev) ! Variables d'etat REAL t(klon, klev) REAL q(klon, klev) ! Profiles cible REAL t_targ(klon, klev) REAL rh_targ(klon, klev) INTEGER k, i REAL zx_qs DO k = 1, klev DO i = 1, klon t_targ(i, k) = t(i, k) IF (t(i, k) 10000.) THEN IF (t(i, k)