--- trunk/libf/phylmd/physiq.f90 2012/07/26 14:37:37 62 +++ trunk/libf/phylmd/physiq.f90 2013/02/18 16:33:12 69 @@ -14,6 +14,7 @@ use aaam_bud_m, only: aaam_bud USE abort_gcm_m, ONLY: abort_gcm + use aeropt_m, only: aeropt use ajsec_m, only: ajsec USE calendar, ONLY: ymds2ju use calltherm_m, only: calltherm @@ -36,6 +37,7 @@ USE dimsoil, ONLY: nsoilmx use drag_noro_m, only: drag_noro USE fcttre, ONLY: foeew, qsatl, qsats, thermcep + use fisrtilp_m, only: fisrtilp USE hgardfou_m, ONLY: hgardfou USE histsync_m, ONLY: histsync USE histwrite_m, ONLY: histwrite @@ -44,6 +46,7 @@ USE ini_histhf_m, ONLY: ini_histhf USE ini_histday_m, ONLY: ini_histday USE ini_histins_m, ONLY: ini_histins + use newmicro_m, only: newmicro USE oasis_m, ONLY: ok_oasis USE orbite_m, ONLY: orbite, zenang USE ozonecm_m, ONLY: ozonecm @@ -53,9 +56,11 @@ USE phytrac_m, ONLY: phytrac USE qcheck_m, ONLY: qcheck use radlwsw_m, only: radlwsw + use readsulfate_m, only: readsulfate use sugwd_m, only: sugwd USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt USE temps, ONLY: annee_ref, day_ref, itau_phy + use unit_nml_m, only: unit_nml USE yoethf_m, ONLY: r2es, rvtmp2 ! Arguments: @@ -116,12 +121,9 @@ logical rnpb parameter(rnpb = .true.) - character(len = 6), save:: ocean + character(len = 6):: ocean = 'force ' ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple") - logical ok_ocean - SAVE ok_ocean - ! "slab" ocean REAL, save:: tslab(klon) ! temperature of ocean slab REAL, save:: seaice(klon) ! glace de mer (kg/m2) @@ -129,13 +131,11 @@ REAL fluxg(klon) ! flux turbulents ocean-atmosphere ! Modele thermique du sol, a activer pour le cycle diurne: - logical, save:: ok_veget - LOGICAL, save:: ok_journe ! sortir le fichier journalier - - LOGICAL ok_mensuel ! sortir le fichier mensuel + logical:: ok_veget = .false. ! type de modele de vegetation utilise - LOGICAL ok_instan ! sortir le fichier instantane - save ok_instan + logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false. + ! sorties journalieres, mensuelles et instantanees dans les + ! fichiers histday, histmth et histins LOGICAL ok_region ! sortir le fichier regional PARAMETER (ok_region = .FALSE.) @@ -318,9 +318,7 @@ SAVE qcondc REAL ema_work1(klon, llm), ema_work2(klon, llm) SAVE ema_work1, ema_work2 - - REAL wd(klon) ! sb - SAVE wd ! sb + REAL, save:: wd(klon) ! Variables locales pour la couche limite (al1): @@ -329,7 +327,7 @@ REAL cdragh(klon) ! drag coefficient pour T and Q REAL cdragm(klon) ! drag coefficient pour vent - !AA Pour phytrac + ! Pour phytrac : REAL ycoefh(klon, llm) ! coef d'echange pour phytrac REAL yu1(klon) ! vents dans la premiere couche U REAL yv1(klon) ! vents dans la premiere couche V @@ -391,7 +389,6 @@ EXTERNAL alboc ! calculer l'albedo sur ocean !KE43 EXTERNAL conema3 ! convect4.3 - EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie) EXTERNAL nuage ! calculer les proprietes radiatives EXTERNAL transp ! transport total de l'eau et de l'energie @@ -449,17 +446,12 @@ REAL dist, rmu0(klon), fract(klon) REAL zdtime ! pas de temps du rayonnement (s) real zlongi - REAL z_avant(klon), z_apres(klon), z_factor(klon) - LOGICAL zx_ajustq - REAL za, zb REAL zx_t, zx_qs, zdelta, zcor real zqsat(klon, llm) INTEGER i, k, iq, nsrf - REAL t_coup - PARAMETER (t_coup = 234.0) - + REAL, PARAMETER:: t_coup = 234. REAL zphi(klon, llm) !IM cf. AM Variables locales pour la CLA (hbtm2) @@ -496,7 +488,6 @@ REAL rflag(klon) ! flag fonctionnement de convect INTEGER iflagctrl(klon) ! flag fonctionnement de convect ! -- convect43: - INTEGER ntra ! nb traceurs pour convect4.3 REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm) REAL dplcldt(klon), dplcldr(klon) @@ -535,20 +526,17 @@ REAL d_u_lif(klon, llm), d_v_lif(klon, llm) REAL d_t_lif(klon, llm) - REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm) - real ratqsbas, ratqshaut - save ratqsbas, ratqshaut, ratqs + REAL, save:: ratqs(klon, llm) + real ratqss(klon, llm), ratqsc(klon, llm) + real:: ratqsbas = 0.01, ratqshaut = 0.3 ! Parametres lies au nouveau schema de nuages (SB, PDF) - real, save:: fact_cldcon - real, save:: facttemps - logical ok_newmicro - save ok_newmicro + real:: fact_cldcon = 0.375 + real:: facttemps = 1.e-4 + logical:: ok_newmicro = .true. real facteur - integer iflag_cldcon - save iflag_cldcon - + integer:: iflag_cldcon = 1 logical ptconv(klon, llm) ! Variables locales pour effectuer les appels en série : @@ -581,8 +569,6 @@ REAL zsto - character(len = 20) modname - character(len = 80) abort_message logical ok_sync real date0 @@ -594,56 +580,52 @@ REAL zero_v(klon) CHARACTER(LEN = 15) tit INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics - INTEGER, SAVE:: if_ebil ! level for energy conservation diagnostics + INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation REAL d_t_ec(klon, llm) ! tendance due à la conversion Ec -> E thermique REAL ZRCPD REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m - REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m - REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille - REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille - !jq Aerosol effects (Johannes Quaas, 27/11/2003) - REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3] + REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m + REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille + REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille + + ! Aerosol effects: + + REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3) REAL, save:: sulfate_pi(klon, llm) - ! (SO4 aerosol concentration, in ug/m3, pre-industrial value) + ! SO4 aerosol concentration, in micro g/m3, pre-industrial value REAL cldtaupi(klon, llm) - ! (Cloud optical thickness for pre-industrial (pi) aerosols) + ! cloud optical thickness for pre-industrial (pi) aerosols REAL re(klon, llm) ! Cloud droplet effective radius REAL fl(klon, llm) ! denominator of re ! Aerosol optical properties - REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2) - REAL cg_ae(klon, llm, 2) - - REAL topswad(klon), solswad(klon) ! Aerosol direct effect. - ! ok_ade = True -ADE = topswad-topsw + REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2) + REAL, save:: cg_ae(klon, llm, 2) + REAL topswad(klon), solswad(klon) ! aerosol direct effect REAL topswai(klon), solswai(klon) ! aerosol indirect effect - ! ok_aie = True -> - ! ok_ade = True -AIE = topswai-topswad - ! ok_ade = F -AIE = topswai-topsw REAL aerindex(klon) ! POLDER aerosol index - ! Parameters - LOGICAL, save:: ok_ade ! apply aerosol direct effect - LOGICAL, save:: ok_aie ! Apply aerosol indirect effect - REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995) + LOGICAL:: ok_ade = .false. ! apply aerosol direct effect + LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect + + REAL:: bl95_b0 = 2., bl95_b1 = 0.2 + ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus + ! B). They link cloud droplet number concentration to aerosol mass + ! concentration. - SAVE bl95_b0, bl95_b1 SAVE u10m SAVE v10m SAVE t2m SAVE q2m SAVE ffonte SAVE fqcalving - SAVE piz_ae - SAVE tau_ae - SAVE cg_ae SAVE rain_con SAVE snow_con SAVE topswai @@ -660,19 +642,17 @@ real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2 + namelist /physiq_nml/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, & + fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, & + ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, & + nsplit_thermals + !---------------------------------------------------------------- - modname = 'physiq' - IF (if_ebil >= 1) THEN - DO i = 1, klon - zero_v(i) = 0. - END DO - END IF + IF (if_ebil >= 1) zero_v = 0. ok_sync = .TRUE. - IF (nqmx < 2) THEN - abort_message = 'eaux vapeur et liquide sont indispensables' - CALL abort_gcm(modname, abort_message, 1) - ENDIF + IF (nqmx < 2) CALL abort_gcm('physiq', & + 'eaux vapeur et liquide sont indispensables', 1) test_firstcal: IF (firstcal) THEN ! initialiser @@ -687,8 +667,6 @@ cg_ae = 0. rain_con(:) = 0. snow_con(:) = 0. - bl95_b0 = 0. - bl95_b1 = 0. topswai(:) = 0. topswad(:) = 0. solswai(:) = 0. @@ -714,11 +692,13 @@ IF (if_ebil >= 1) d_h_vcol_phy = 0. - ! Appel à la lecture du run.def physique - call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, & - fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, & - ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, & - iflag_thermals, nsplit_thermals) + iflag_thermals = 0 + nsplit_thermals = 1 + print *, "Enter namelist 'physiq_nml'." + read(unit=*, nml=physiq_nml) + write(unit_nml, nml=physiq_nml) + + call conf_phys ! Initialiser les compteurs: @@ -732,7 +712,7 @@ ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0) ! ATTENTION : il faudra a terme relire q2 dans l'etat initial - q2 = 1.e-8 + q2 = 1e-8 radpas = NINT(86400. / dtphys / nbapp_rad) @@ -740,30 +720,19 @@ IF (raz_date) itau_phy = 0 PRINT *, 'cycle_diurne = ', cycle_diurne + CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, & + ok_instan, ok_region) - IF(ocean.NE.'force ') THEN - ok_ocean = .TRUE. - ENDIF - - CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, & - ok_region) - - IF (dtphys*REAL(radpas) > 21600..AND.cycle_diurne) THEN - print *, 'Nbre d appels au rayonnement insuffisant' + IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN print *, "Au minimum 4 appels par jour si cycle diurne" - abort_message = 'Nbre d appels au rayonnement insuffisant' - call abort_gcm(modname, abort_message, 1) + call abort_gcm('physiq', & + "Nombre d'appels au rayonnement insuffisant", 1) ENDIF - print *, "Clef pour la convection, iflag_con = ", iflag_con - ! Initialisation pour la convection de K.E. (sb): + ! Initialisation pour le schéma de convection d'Emanuel : IF (iflag_con >= 3) THEN - print *, "Convection de Kerry Emanuel 4.3" - - DO i = 1, klon - ibas_con(i) = 1 - itop_con(i) = 1 - ENDDO + ibas_con = 1 + itop_con = 1 ENDIF IF (ok_orodr) THEN @@ -793,19 +762,19 @@ call ini_histday(dtphys, ok_journe, nid_day, nqmx) call ini_histins(dtphys, ok_instan, nid_ins) CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0) - !XXXPB Positionner date0 pour initialisation de ORCHIDEE - WRITE(*, *) 'physiq date0: ', date0 + ! Positionner date0 pour initialisation de ORCHIDEE + print *, 'physiq date0: ', date0 ENDIF test_firstcal ! Mettre a zero des variables de sortie (pour securite) DO i = 1, klon - d_ps(i) = 0.0 + d_ps(i) = 0. ENDDO DO iq = 1, nqmx DO k = 1, llm DO i = 1, klon - d_qx(i, k, iq) = 0.0 + d_qx(i, k, iq) = 0. ENDDO ENDDO ENDDO @@ -890,7 +859,7 @@ forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg - ! Mettre en action les conditions aux limites (albedo, sst, etc.). + ! Mettre en action les conditions aux limites (albedo, sst etc.). ! Prescrire l'ozone et calculer l'albedo sur l'ocean. wo = ozonecm(REAL(julien), paprs) @@ -959,8 +928,8 @@ DO nsrf = 1, nbsrf DO i = 1, klon fsollw(i, nsrf) = sollw(i) & - + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf)) - fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i)) + + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf)) + fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i)) ENDDO ENDDO @@ -1001,7 +970,7 @@ END DO DO i = 1, klon sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol - evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol + evap(i) = - zxfluxq(i, 1) ! flux d'évaporation au sol fder(i) = dlw(i) + dsens(i) + devap(i) ENDDO @@ -1048,12 +1017,9 @@ s_trmb2(i) = 0.0 s_trmb3(i) = 0.0 - IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + & - pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) > EPSFRA) & - THEN - WRITE(*, *) 'physiq : pb sous surface au point ', i, & - pctsrf(i, 1 : nbsrf) - ENDIF + IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) & + + pctsrf(i, is_sic) - 1.) > EPSFRA) print *, & + 'physiq : problème sous surface au point ', i, pctsrf(i, 1 : nbsrf) ENDDO DO nsrf = 1, nbsrf DO i = 1, klon @@ -1110,38 +1076,25 @@ ! Calculer la derive du flux infrarouge DO i = 1, klon - dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3 + dlw(i) = - 4. * RSIGMA * zxtsol(i)**3 ENDDO ! Appeler la convection (au choix) DO k = 1, llm DO i = 1, klon - conv_q(i, k) = d_q_dyn(i, k) & - + d_q_vdf(i, k)/dtphys - conv_t(i, k) = d_t_dyn(i, k) & - + d_t_vdf(i, k)/dtphys + conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys + conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k)/dtphys ENDDO ENDDO + IF (check) THEN za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy) print *, "avantcon = ", za ENDIF - zx_ajustq = iflag_con == 2 - IF (zx_ajustq) THEN - DO i = 1, klon - z_avant(i) = 0.0 - ENDDO - DO k = 1, llm - DO i = 1, klon - z_avant(i) = z_avant(i) + (q_seri(i, k) + ql_seri(i, k)) & - *zmasse(i, k) - ENDDO - ENDDO - ENDIF - select case (iflag_con) - case (2) + if (iflag_con == 2) then + z_avant = sum((q_seri + ql_seri) * zmasse, dim=2) CALL conflx(dtphys, paprs, play, t_seri, q_seri, conv_t, conv_q, & zxfluxq(1, 1), omega, d_t_con, d_q_con, rain_con, snow_con, pmfu, & pmfd, pen_u, pde_u, pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, & @@ -1152,29 +1105,22 @@ ibas_con(i) = llm + 1 - kcbot(i) itop_con(i) = llm + 1 - kctop(i) ENDDO - case (3:) - ! number of tracers for the convection scheme of Kerry Emanuel: + else + ! iflag_con >= 3 + CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, & + v_seri, tr_seri, ema_work1, ema_work2, d_t_con, d_q_con, & + d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, & + itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, & + pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, & + wd, pmflxr, pmflxs, da, phi, mp, ntra=1) + ! (number of tracers for the convection scheme of Kerry Emanuel: ! la partie traceurs est faite dans phytrac ! on met ntra = 1 pour limiter les appels mais on peut - ! supprimer les calculs / ftra. - ntra = 1 - ! Schéma de convection modularisé et vectorisé : - ! (driver commun aux versions 3 et 4) - - CALL concvl(iflag_con, dtphys, paprs, play, t_seri, q_seri, u_seri, & - v_seri, tr_seri, ntra, ema_work1, ema_work2, d_t_con, d_q_con, & - d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, itop_con, & - upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, pbase, bbase, & - dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, wd, pmflxr, pmflxs, & - da, phi, mp) + ! supprimer les calculs / ftra.) + clwcon0 = qcondc pmfu = upwd + dnwd - - IF (.NOT. ok_gust) THEN - do i = 1, klon - wd(i) = 0.0 - enddo - ENDIF + IF (.NOT. ok_gust) wd = 0. ! Calcul des propriétés des nuages convectifs @@ -1202,10 +1148,7 @@ clwcon0 = fact_cldcon*clwcon0 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, & rnebcon0) - case default - print *, "iflag_con non-prevu", iflag_con - stop 1 - END select + END if DO k = 1, llm DO i = 1, klon @@ -1239,20 +1182,10 @@ zx_t = zx_t/za*dtphys print *, "Precip = ", zx_t ENDIF - IF (zx_ajustq) THEN - DO i = 1, klon - z_apres(i) = 0.0 - ENDDO - DO k = 1, llm - DO i = 1, klon - z_apres(i) = z_apres(i) + (q_seri(i, k) + ql_seri(i, k)) & - *zmasse(i, k) - ENDDO - ENDDO - DO i = 1, klon - z_factor(i) = (z_avant(i)-(rain_con(i) + snow_con(i))*dtphys) & - /z_apres(i) - ENDDO + + IF (iflag_con == 2) THEN + z_apres = sum((q_seri + ql_seri) * zmasse, dim=2) + z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres DO k = 1, llm DO i = 1, klon IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN @@ -1261,7 +1194,6 @@ ENDDO ENDDO ENDIF - zx_ajustq = .FALSE. ! Convection sèche (thermiques ou ajustement) @@ -1316,7 +1248,7 @@ enddo ! ratqs final - if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then + if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then ! les ratqs sont une conbinaison de ratqss et ratqsc ! ratqs final ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de @@ -1393,9 +1325,8 @@ endif ! Nuages diagnostiques pour Tiedtke - CALL diagcld1(paprs, play, & - rain_tiedtke, snow_tiedtke, ibas_con, itop_con, & - diafra, dialiq) + CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, & + itop_con, diafra, dialiq) DO k = 1, llm DO i = 1, klon IF (diafra(i, k) > cldfra(i, k)) THEN @@ -1472,13 +1403,11 @@ ENDDO ! Introduce the aerosol direct and first indirect radiative forcings: - ! Johannes Quaas, 27/11/2003 IF (ok_ade .OR. ok_aie) THEN - ! Get sulfate aerosol distribution + ! Get sulfate aerosol distribution : CALL readsulfate(rdayvrai, firstcal, sulfate) CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi) - ! Calculate aerosol optical properties (Olivier Boucher) CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, & aerindex) ELSE @@ -1489,10 +1418,9 @@ ! Paramètres optiques des nuages et quelques paramètres pour diagnostics : if (ok_newmicro) then - CALL newmicro(paprs, play, ok_newmicro, t_seri, cldliq, cldfra, & - cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, & - fiwc, ok_aie, sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, & - re, fl) + CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, & + cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, & + sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl) else CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, & cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, & @@ -1688,7 +1616,7 @@ ! SORTIES - !cc prw = eau precipitable + ! prw = eau precipitable DO i = 1, klon prw(i) = 0. DO k = 1, llm