--- trunk/libf/phylmd/clmain.f90 2011/01/06 17:52:19 38 +++ trunk/libf/phylmd/clmain.f90 2011/07/01 15:00:48 47 @@ -16,22 +16,22 @@ fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice) ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19 + ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18 + ! Objet : interface de "couche limite" (diffusion verticale) - ! Tout ce qui a trait aux traceurs est dans phytrac maintenant. + ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant. ! Pour l'instant le calcul de la couche limite pour les traceurs - ! se fait avec cltrac et ne tient pas compte de la différentiation + ! se fait avec "cltrac" et ne tient pas compte de la différentiation ! des sous-fractions de sol. ! Pour pouvoir extraire les coefficients d'échanges et le vent - ! dans la première couche, trois champs supplémentaires ont été créés : - ! zcoefh, zu1 et zv1. Pour l'instant nous avons moyenné les valeurs - ! de ces trois champs sur les 4 sous-surfaces du modèle. Dans l'avenir - ! si les informations des sous-surfaces doivent être prises en compte - ! il faudra sortir ces mêmes champs en leur ajoutant une dimension, - ! c'est a dire nbsrf (nombre de sous-surfaces). - - ! Auteur Z.X. Li (LMD/CNRS) date: 1993/08/18 - ! Objet : interface de "couche limite" (diffusion verticale) + ! dans la première couche, trois champs supplémentaires ont été + ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons + ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces + ! du modèle. Dans l'avenir, si les informations des sous-surfaces + ! doivent être prises en compte, il faudra sortir ces mêmes champs + ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de + ! sous-surfaces). ! Arguments: ! dtime----input-R- interval du temps (secondes) @@ -92,29 +92,32 @@ ! pblh------- HCL ! pblT------- T au nveau HCL - USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync - use histwrite_m, only: histwrite use calendar, ONLY : ymds2ju + use coefkz_m, only: coefkz + use coefkzmin_m, only: coefkzmin + USE conf_phys_m, ONLY : iflag_pbl USE dimens_m, ONLY : iim, jjm - USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf USE dimphy, ONLY : klev, klon, zmasq USE dimsoil, ONLY : nsoilmx - USE temps, ONLY : annee_ref, itau_phy USE dynetat0_m, ONLY : day_ini - USE iniprint, ONLY : prt_level - USE suphec_m, ONLY : rd, rg, rkappa - USE conf_phys_m, ONLY : iflag_pbl USE gath_cpl, ONLY : gath2cpl use hbtm_m, only: hbtm + USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync + use histwrite_m, only: histwrite + USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf + USE iniprint, ONLY : prt_level + USE suphec_m, ONLY : rd, rg, rkappa + USE temps, ONLY : annee_ref, itau_phy + use yamada4_m, only: yamada4 REAL, INTENT (IN) :: dtime REAL date0 INTEGER, INTENT (IN) :: itap REAL t(klon, klev), q(klon, klev) - REAL u(klon, klev), v(klon, klev) - REAL, INTENT (IN) :: paprs(klon, klev+1) - REAL, INTENT (IN) :: pplay(klon, klev) - REAL, INTENT (IN) :: rlon(klon), rlat(klon) + REAL, INTENT (IN):: u(klon, klev), v(klon, klev) + REAL, INTENT (IN):: paprs(klon, klev+1) + REAL, INTENT (IN):: pplay(klon, klev) + REAL, INTENT (IN):: rlon(klon), rlat(klon) REAL cufi(klon), cvfi(klon) REAL d_t(klon, klev), d_q(klon, klev) REAL d_u(klon, klev), d_v(klon, klev) @@ -179,7 +182,7 @@ REAL ytsoil(klon, nsoilmx) REAL qsol(klon) - EXTERNAL clqh, clvent, coefkz, calbeta, cltrac + EXTERNAL clqh, clvent, calbeta, cltrac REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon) REAL yalb(klon) @@ -208,22 +211,21 @@ PARAMETER (ok_nonloc=.FALSE.) REAL ycoefm0(klon, klev), ycoefh0(klon, klev) - !IM 081204 hcl_Anne ? BEG REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev) REAL ykmm(klon, klev+1), ykmn(klon, klev+1) REAL ykmq(klon, klev+1) REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf) REAL q2diag(klon, klev+1) - !IM 081204 hcl_Anne ? END REAL u1lay(klon), v1lay(klon) REAL delp(klon, klev) INTEGER i, k, nsrf INTEGER ni(klon), knon, j - ! Introduction d'une variable "pourcentage potentiel" pour tenir compte - ! des eventuelles apparitions et/ou disparitions de la glace de mer + REAL pctsrf_pot(klon, nbsrf) + ! "pourcentage potentiel" pour tenir compte des éventuelles + ! apparitions ou disparitions de la glace de mer REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola. @@ -298,7 +300,6 @@ !------------------------------------------------------------ - ! initialisation Anne ytherm = 0. IF (debugindex .AND. first_appel) THEN @@ -307,7 +308,7 @@ ! initialisation sorties netcdf idayref = day_ini - CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian) + CALL ymds2ju(annee_ref, 1, idayref, 0., zjulian) CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon) DO i = 1, iim zx_lon(i, 1) = rlon(i+1) @@ -342,75 +343,64 @@ v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2 END DO - ! initialisation: - - DO i = 1, klon - rugmer(i) = 0.0 - cdragh(i) = 0.0 - cdragm(i) = 0.0 - dflux_t(i) = 0.0 - dflux_q(i) = 0.0 - zu1(i) = 0.0 - zv1(i) = 0.0 - END DO - ypct = 0.0 - yts = 0.0 - ysnow = 0.0 - yqsurf = 0.0 - yalb = 0.0 - yalblw = 0.0 - yrain_f = 0.0 - ysnow_f = 0.0 - yfder = 0.0 - ytaux = 0.0 - ytauy = 0.0 - ysolsw = 0.0 - ysollw = 0.0 - ysollwdown = 0.0 - yrugos = 0.0 - yu1 = 0.0 - yv1 = 0.0 - yrads = 0.0 - ypaprs = 0.0 - ypplay = 0.0 - ydelp = 0.0 - yu = 0.0 - yv = 0.0 - yt = 0.0 - yq = 0.0 - pctsrf_new = 0.0 - y_flux_u = 0.0 - y_flux_v = 0.0 + ! Initialization: + rugmer = 0. + cdragh = 0. + cdragm = 0. + dflux_t = 0. + dflux_q = 0. + zu1 = 0. + zv1 = 0. + ypct = 0. + yts = 0. + ysnow = 0. + yqsurf = 0. + yalb = 0. + yalblw = 0. + yrain_f = 0. + ysnow_f = 0. + yfder = 0. + ytaux = 0. + ytauy = 0. + ysolsw = 0. + ysollw = 0. + ysollwdown = 0. + yrugos = 0. + yu1 = 0. + yv1 = 0. + yrads = 0. + ypaprs = 0. + ypplay = 0. + ydelp = 0. + yu = 0. + yv = 0. + yt = 0. + yq = 0. + pctsrf_new = 0. + y_flux_u = 0. + y_flux_v = 0. !$$ PB - y_dflux_t = 0.0 - y_dflux_q = 0.0 + y_dflux_t = 0. + y_dflux_q = 0. ytsoil = 999999. yrugoro = 0. ! -- LOOP - yu10mx = 0.0 - yu10my = 0.0 - ywindsp = 0.0 + yu10mx = 0. + yu10my = 0. + ywindsp = 0. ! -- LOOP - DO nsrf = 1, nbsrf - DO i = 1, klon - d_ts(i, nsrf) = 0.0 - END DO - END DO + d_ts = 0. !§§§ PB yfluxlat = 0. flux_t = 0. flux_q = 0. flux_u = 0. flux_v = 0. - DO k = 1, klev - DO i = 1, klon - d_t(i, k) = 0.0 - d_q(i, k) = 0.0 - d_u(i, k) = 0.0 - d_v(i, k) = 0.0 - zcoefh(i, k) = 0.0 - END DO - END DO + d_t = 0. + d_q = 0. + d_u = 0. + d_v = 0. + zcoefh = 0. ! Boucler sur toutes les sous-fractions du sol: @@ -427,7 +417,7 @@ ni = 0 knon = 0 DO i = 1, klon - ! pour determiner le domaine a traiter on utilise les surfaces + ! Pour déterminer le domaine à traiter, on utilise les surfaces ! "potentielles" IF (pctsrf_pot(i, nsrf) > epsfra) THEN knon = knon + 1 @@ -447,7 +437,7 @@ CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab) END IF - IF (knon==0) CYCLE + IF (knon == 0) CYCLE DO j = 1, knon i = ni(j) @@ -479,8 +469,8 @@ ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j)) END DO - ! IF bucket model for continent, copy soil water content - IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN + ! IF bucket model for continent, copy soil water content + IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN DO j = 1, knon i = ni(j) yqsol(j) = qsol(i) @@ -511,10 +501,7 @@ ! calculer Cdrag et les coefficients d'echange CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,& yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh) - !IM 081204 BEG - !CR test - IF (iflag_pbl==1) THEN - !IM 081204 END + IF (iflag_pbl == 1) THEN CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0) DO k = 1, klev DO i = 1, knon @@ -524,36 +511,28 @@ END DO END IF - !IM cf JLD : on seuille ycoefm et ycoefh - IF (nsrf==is_oce) THEN + ! on seuille ycoefm et ycoefh + IF (nsrf == is_oce) THEN DO j = 1, knon - ! ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3) ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax) - ! ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3) ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax) END DO END IF - !IM: 261103 IF (ok_kzmin) THEN - !IM cf FH: 201103 BEG - ! Calcul d'une diffusion minimale pour les conditions tres stables. - CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, & + ! Calcul d'une diffusion minimale pour les conditions tres stables + CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm(:, 1), & ycoefm0, ycoefh0) - IF (1==1) THEN - DO k = 1, klev - DO i = 1, knon - ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k)) - ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k)) - END DO + DO k = 1, klev + DO i = 1, knon + ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k)) + ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k)) END DO - END IF - !IM cf FH: 201103 END - !IM: 261103 - END IF !ok_kzmin + END DO + END IF - IF (iflag_pbl>=3) THEN + IF (iflag_pbl >= 3) THEN ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, & 1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg @@ -579,30 +558,23 @@ END DO END DO - ! Bug introduit volontairement pour converger avec les resultats - ! du papier sur les thermiques. - IF (1==1) THEN - y_cd_m(1:knon) = ycoefm(1:knon, 1) - y_cd_h(1:knon) = ycoefh(1:knon, 1) - ELSE - y_cd_h(1:knon) = ycoefm(1:knon, 1) - y_cd_m(1:knon) = ycoefh(1:knon, 1) - END IF + y_cd_m(1:knon) = ycoefm(1:knon, 1) + y_cd_h(1:knon) = ycoefh(1:knon, 1) CALL ustarhb(knon, yu, yv, y_cd_m, yustar) IF (prt_level>9) THEN PRINT *, 'USTAR = ', yustar END IF - ! iflag_pbl peut etre utilise comme longuer de melange + ! iflag_pbl peut être utilisé comme longueur de mélange - IF (iflag_pbl>=11) THEN + IF (iflag_pbl >= 11) THEN CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, & yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, & iflag_pbl) ELSE - CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, & - yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl) + CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, & + y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl) END IF ycoefm(1:knon, 1) = y_cd_m(1:knon) @@ -621,12 +593,6 @@ ytaux = y_flux_u(:, 1) ytauy = y_flux_v(:, 1) - ! FH modif sur le cdrag temperature - !$$$PB : déplace dans clcdrag - !$$$ do i=1, knon - !$$$ ycoefh(i, 1)=ycoefm(i, 1)*0.8 - !$$$ enddo - ! calculer la diffusion de "q" et de "h" CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,& cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,& @@ -641,7 +607,7 @@ ! calculer la longueur de rugosite sur ocean yrugm = 0. - IF (nsrf==is_oce) THEN + IF (nsrf == is_oce) THEN DO j = 1, knon yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + & 0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)) @@ -662,17 +628,12 @@ ycoefm(j, k) = ycoefm(j, k)*ypct(j) y_d_t(j, k) = y_d_t(j, k)*ypct(j) y_d_q(j, k) = y_d_q(j, k)*ypct(j) - !§§§ PB flux_t(i, k, nsrf) = y_flux_t(j, k) flux_q(i, k, nsrf) = y_flux_q(j, k) flux_u(i, k, nsrf) = y_flux_u(j, k) flux_v(i, k, nsrf) = y_flux_v(j, k) - !$$$ PB y_flux_t(j, k) = y_flux_t(j, k) * ypct(j) - !$$$ PB y_flux_q(j, k) = y_flux_q(j, k) * ypct(j) y_d_u(j, k) = y_d_u(j, k)*ypct(j) y_d_v(j, k) = y_d_v(j, k)*ypct(j) - !$$$ PB y_flux_u(j, k) = y_flux_u(j, k) * ypct(j) - !$$$ PB y_flux_v(j, k) = y_flux_v(j, k) * ypct(j) END DO END DO @@ -693,12 +654,10 @@ qsurf(i, nsrf) = yqsurf(j) rugos(i, nsrf) = yz0_new(j) fluxlat(i, nsrf) = yfluxlat(j) - !$$$ pb rugmer(i) = yrugm(j) - IF (nsrf==is_oce) THEN + IF (nsrf == is_oce) THEN rugmer(i) = yrugm(j) rugos(i, nsrf) = yrugm(j) END IF - !IM cf JLD ?? agesno(i, nsrf) = yagesno(j) fqcalving(i, nsrf) = y_fqcalving(j) ffonte(i, nsrf) = y_ffonte(j) @@ -709,13 +668,13 @@ zu1(i) = zu1(i) + yu1(j) zv1(i) = zv1(i) + yv1(j) END DO - IF (nsrf==is_ter) THEN + IF (nsrf == is_ter) THEN DO j = 1, knon i = ni(j) qsol(i) = yqsol(j) END DO END IF - IF (nsrf==is_lic) THEN + IF (nsrf == is_lic) THEN DO j = 1, knon i = ni(j) run_off_lic_0(i) = y_run_off_lic_0(j) @@ -735,12 +694,8 @@ DO k = 1, klev d_t(i, k) = d_t(i, k) + y_d_t(j, k) d_q(i, k) = d_q(i, k) + y_d_q(j, k) - !$$$ PB flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k) - !$$$ flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k) d_u(i, k) = d_u(i, k) + y_d_u(j, k) d_v(i, k) = d_v(i, k) + y_d_v(j, k) - !$$$ PB flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k) - !$$$ flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k) zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k) END DO END DO @@ -757,7 +712,7 @@ 1)))*(ypaprs(j, 1)-ypplay(j, 1)) tairsol(j) = yts(j) + y_d_ts(j) rugo1(j) = yrugos(j) - IF (nsrf==is_oce) THEN + IF (nsrf == is_oce) THEN rugo1(j) = rugos(i, nsrf) END IF psfce(j) = ypaprs(j, 1) @@ -769,7 +724,6 @@ CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, & tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, & yu10m, yustar) - !IM 081204 END DO j = 1, knon i = ni(j) @@ -811,7 +765,7 @@ END DO END DO !IM "slab" ocean - IF (nsrf==is_oce) THEN + IF (nsrf == is_oce) THEN DO j = 1, knon ! on projette sur la grille globale i = ni(j) @@ -823,11 +777,10 @@ END DO END IF - IF (nsrf==is_sic) THEN + IF (nsrf == is_sic) THEN DO j = 1, knon i = ni(j) ! On pondère lorsque l'on fait le bilan au sol : - ! flux_g(i) = y_flux_g(j)*ypct(j) IF (pctsrf_new(i, is_sic)>epsfra) THEN flux_g(i) = y_flux_g(j) ELSE @@ -836,19 +789,15 @@ END DO END IF - !nsrf.EQ.is_sic - IF (ocean=='slab ') THEN - IF (nsrf==is_oce) THEN + IF (ocean == 'slab ') THEN + IF (nsrf == is_oce) THEN tslab(1:klon) = ytslab(1:klon) seaice(1:klon) = y_seaice(1:klon) - !nsrf END IF - !OCEAN END IF END DO ! On utilise les nouvelles surfaces - ! A rajouter: conservation de l'albedo rugos(:, is_oce) = rugmer pctsrf = pctsrf_new