Index: LMDZ5/trunk/libf/dyn3d/gcm.F =================================================================== --- LMDZ5/trunk/libf/dyn3d/gcm.F (revision 1527) +++ LMDZ5/trunk/libf/dyn3d/gcm.F (revision 1529) @@ -408,5 +408,5 @@ c ------------------------------- - IF (call_iniphys.and.(iflag_phys.eq.1)) THEN + IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN latfi(1)=rlatu(1) lonfi(1)=0. @@ -459,6 +459,15 @@ #endif +#ifdef CPP_EARTH +! Create start file (startphy.nc) and boundary conditions (limit.nc) +! for the Earth verstion + if (iflag_phys>=100) then + call iniaqua(ngridmx,latfi,lonfi,iflag_phys) + endif +#endif + ! if (planet_type.eq."earth") then ! Write an Earth-format restart file + CALL dynredem0("restart.nc", day_end, phis) ! endif Index: LMDZ5/trunk/libf/dyn3d/iniacademic.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3d/iniacademic.F90 (revision 1527) +++ LMDZ5/trunk/libf/dyn3d/iniacademic.F90 (revision 1529) @@ -115,5 +115,5 @@ endif - academic_case: if (iflag_phys == 2) then + academic_case: if (iflag_phys >= 2) then ! initializations @@ -208,5 +208,9 @@ IF (.NOT. read_start) THEN ! surface pressure - ps(:)=preff + if (iflag_phys>2) then + ps(:)=preff + else + ps(:)=101080. + endif ! ground geopotential phis(:)=0. Index: LMDZ5/trunk/libf/dyn3d/leapfrog.F =================================================================== --- LMDZ5/trunk/libf/dyn3d/leapfrog.F (revision 1527) +++ LMDZ5/trunk/libf/dyn3d/leapfrog.F (revision 1529) @@ -149,4 +149,5 @@ logical ok_sync parameter (ok_sync = .true.) + logical physic data callinigrads/.true./ @@ -193,7 +194,8 @@ itaufin = nday*day_step itaufinp1 = itaufin +1 - - itau = 0 + physic=.true. + if (iflag_phys==0.or.iflag_phys==2) physic=.false. + c iday = day_ini+itau/day_step c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0 @@ -281,5 +283,5 @@ s apdiss = .TRUE. IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward - s .and. iflag_phys.EQ.1 ) apphys = .TRUE. + s .and. physic ) apphys = .TRUE. ELSE ! Leapfrog/Matsuno time stepping @@ -287,5 +289,5 @@ IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward ) s apdiss = .TRUE. - IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE. + IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE. END IF Index: LMDZ5/trunk/libf/phylmd/phyaqua.F =================================================================== --- LMDZ5/trunk/libf/phylmd/phyaqua.F (revision 1529) +++ LMDZ5/trunk/libf/phylmd/phyaqua.F (revision 1529) @@ -0,0 +1,757 @@ +! Routines complementaires pour la physique planetaire. + + + subroutine iniaqua(nlon,latfi,lonfi,iflag_phys) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Creation d'un etat initial et de conditions aux limites +! (resp startphy.nc et limit.nc) pour des configurations idealisees +! du modele LMDZ dans sa version terrestre. +! iflag_phys est un parametre qui controle +! iflag_phys = N +! de 100 a 199 : aqua planetes avec SST forcees +! N-100 determine le type de SSTs +! de 200 a 299 : terra planetes avec Ts calcule +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + use comgeomphy + use dimphy + use surface_data, only : type_ocean,ok_veget + use pbl_surface_mod, only : pbl_surface_init + USE fonte_neige_mod, only : fonte_neige_init + use phys_state_var_mod + use control_mod + + + USE IOIPSL + IMPLICIT NONE + +#include "dimensions.h" +! #include "dimphy.h" +! #include "YOMCST.h" +#include "comconst.h" +#include "clesphys.h" +#include "dimsoil.h" +#include "indicesol.h" + + integer nlon,iflag_phys +cIM ajout latfi, lonfi + REAL, DIMENSION (nlon) :: lonfi, latfi + INTEGER type_profil,type_aqua + +c Ajouts initialisation des surfaces + REAL :: run_off_lic_0(nlon) + REAL :: qsolsrf(nlon,nbsrf),snsrf(nlon,nbsrf) + REAL :: frugs(nlon,nbsrf) + REAL :: agesno(nlon,nbsrf) + REAL :: tsoil(nlon,nsoilmx,nbsrf) + REAL :: tslab(nlon), seaice(nlon) + REAL evap(nlon,nbsrf),fder(nlon) + + + +c Arguments : +c ----------- + +! integer radpas + integer it,unit,i,k,itap + + real airefi,zcufi,zcvfi + + real rugos,albedo + REAL tsurf + REAL time,timestep,day,day0 + real qsol_f,qsol(nlon) + real rugsrel(nlon) +! real zmea(nlon),zstd(nlon),zsig(nlon) +! real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon) +! real rlon(nlon),rlat(nlon) + logical alb_ocean +! integer demih_pas + + integer day_ini + + CHARACTER*80 ans,file_forctl, file_fordat, file_start + character*100 file,var + character*2 cnbl + + REAL phy_nat(nlon,360) + REAL phy_alb(nlon,360) + REAL phy_sst(nlon,360) + REAL phy_bil(nlon,360) + REAL phy_rug(nlon,360) + REAL phy_ice(nlon,360) + REAL phy_fter(nlon,360) + REAL phy_foce(nlon,360) + REAL phy_fsic(nlon,360) + REAL phy_flic(nlon,360) + + integer, save:: read_climoz ! read ozone climatology + + +!------------------------------------------------------------------------- +! declaration pour l'appel a phyredem +!------------------------------------------------------------------------- + +! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf) + real falbe(nlon,nbsrf),falblw(nlon,nbsrf) +! real pbl_tke(nlon,llm,nbsrf) +! real rain_fall(nlon),snow_fall(nlon) +! real solsw(nlon), sollw(nlon),radsol(nlon) +! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm) +! real ratqs(nlon,llm) +! real clwcon(nlon,llm) + + INTEGER longcles + PARAMETER ( longcles = 20 ) + REAL clesphy0( longcles ) + + +c----------------------------------------------------------------------- +c dynamial tendencies : +c --------------------- + + INTEGER l,ierr,aslun + + REAL longitude,latitude + REAL paire + + DATA latitude,longitude/48.,0./ + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! INITIALISATIONS +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!----------------------------------------------------------------------- +! Initialisations des constantes +! ------------------------------- + + + type_aqua=iflag_phys/100 + type_profil=iflag_phys-type_aqua*100 + print*,'type_aqua, type_profil',type_aqua, type_profil + + if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua' + call phys_state_var_init(read_climoz) + + + read_climoz=0 + day0=217. + day=day0 + it=0 + time=0. + +cIM ajout latfi, lonfi + rlatd=latfi + rlond=lonfi + rlat=rlatd*180./pi + rlon=rlond*180./pi + +!----------------------------------------------------------------------- +! initialisations de la physique +!----------------------------------------------------------------------- + + day_ini=dayref + airefi=1. + zcufi=1. + zcvfi=1. + nbapp_rad=24 + CALL getin('nbapp_rad',nbapp_rad) + +!--------------------------------------------------------------------- +c Creation des conditions aux limites: +c ------------------------------------ +! Initialisations des constantes +! Ajouter les manquants dans planete.def... (albedo etc) + co2_ppm=348. + CALL getin('co2_ppm',co2_ppm) + solaire=1365. + CALL getin('solaire',solaire) + radsol=0. + qsol_f=10. + CALL getin('albedo',albedo) + alb_ocean=.true. + CALL getin('alb_ocean',alb_ocean) + +c Conditions aux limites: +c ----------------------- + + qsol(:) = qsol_f + rugsrel = 0.0 ! (rugsrel = rugoro) + agesno = 50.0 +! Relief plat + zmea = 0. + zstd = 0. + zsig = 0. + zgam = 0. + zthe = 0. + zpic = 0. + zval = 0. + +! Une seule surface + pctsrf=0. + if (type_aqua==1) then + rugos=1.e-4 + albedo=0.19 + pctsrf(:,is_oce)=1. + else if (type_aqua==2) then + rugos=0.03 + albedo=0.1 + pctsrf(:,is_ter)=1. + endif + + CALL getin('rugos',rugos) + zmasq(:)=pctsrf(:,is_oce) + +! pctsrf_pot(:,is_oce) = 1. - zmasq(:) +! pctsrf_pot(:,is_sic) = 1. - zmasq(:) + +! Si alb_ocean on calcule un albedo oceanique moyen +! if (alb_ocean) then +! Voir pourquoi on avait ca. +! CALL ini_alb_oce(phy_alb) +! else + phy_alb(:,:) = albedo ! albedo land only (old value condsurf_jyg=0.3) +! endif !alb_ocean + + do i=1,360 +cIM Terraplanete phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2 +cIM ajout calcul profil sst selon le cas considere (cf. FBr) + + phy_nat(:,i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise + phy_bil(:,i) = 1.0 ! ne sert que pour les slab_ocean + phy_rug(:,i) = rugos ! longueur rugosite utilisee sur land only + phy_ice(:,i) = 0.0 ! fraction de glace (?) + phy_fter(:,i) = pctsrf(:,is_ter) ! fraction de glace (?) + phy_foce(:,i) = pctsrf(:,is_oce) ! fraction de glace (?) + phy_fsic(:,i) = pctsrf(:,is_sic) ! fraction de glace (?) + phy_flic(:,i) = pctsrf(:,is_lic) ! fraction de glace (?) + enddo +cIM calcul profil sst + call profil_sst(nlon, rlatd, type_profil, phy_sst) + + call writelim + s (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice, + s phy_fter,phy_foce,phy_flic,phy_fsic) + + +!--------------------------------------------------------------------- +c Ecriture de l'etat initial: +c --------------------------- + +C +C Ecriture etat initial physique +C + timestep = dtvr * FLOAT(iphysiq) + radpas = NINT (daysec/timestep/ FLOAT(nbapp_rad) ) + + DO i = 1, longcles + clesphy0(i) = 0. + ENDDO + clesphy0(1) = FLOAT( iflag_con ) + clesphy0(2) = FLOAT( nbapp_rad ) +c IF( cycle_diurne ) clesphy0(3) = 1. + clesphy0(3)=1. ! cycle_diurne + clesphy0(4)=1. ! soil_model + clesphy0(5)=1. ! new_oliq + clesphy0(6)=0. ! ok_orodr + clesphy0(7)=0. ! ok_orolf + clesphy0(8)=0. ! ok_limitvrai + + +c======================================================================= +c Profils initiaux +c======================================================================= + +! On initialise les temperatures de surfaces comme les sst + do i=1,nlon + ftsol(i,:)=phy_sst(i,1) + tsoil(i,:,:)=phy_sst(i,1) + tslab(i)=phy_sst(i,1) + enddo + + falbe(:,:)=albedo + falblw(:,:)=albedo + rain_fall(:)=0. + snow_fall(:)=0. + solsw(:)=0. + sollw(:)=0. + radsol(:)=0. + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! intialisation bidon mais pas grave + t_ancien(:,:)=0. + q_ancien(:,:)=0. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + rnebcon=0. + ratqs=0. + clwcon=0. + pbl_tke=1.e-8 + +! variables supplementaires pour appel a plb_surface_init + fder(:)=0. + seaice(:)=0. + run_off_lic_0=0. + evap=0. + + +! Initialisations necessaires avant phyredem + type_ocean = "force" + call fonte_neige_init(run_off_lic_0) + qsolsrf(:,:)=qsol(1) ! humidite du sol des sous surface + snsrf(:,:)=0. ! couverture de neige des sous surface + frugs(:,:)=rugos ! couverture de neige des sous surface + + + call pbl_surface_init(qsol, fder, snsrf, qsolsrf, + . evap, frugs, agesno, tsoil) + + print*,'avant phyredem dans iniaqua' + + CALL phyredem ("startphy.nc") + + print*,'apres phyredem' + call phys_state_var_end + + return + end + + +c==================================================================== +c==================================================================== + SUBROUTINE zenang_an(cycle_diurne,gmtime,rlat,rlon,rmu0,fract) + USE dimphy + IMPLICIT none +c==================================================================== +c============================================================= +c CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract) +c Auteur : A. Campoy et F. Hourdin +c Objet : calculer les valeurs moyennes du cos de l'angle zenithal +c et l'ensoleillement moyen entre gmtime1 et gmtime2 +c connaissant la declinaison, la latitude et la longitude. +c +c Dans cette version particuliere, on calcule le rayonnement +c moyen sur l'année à chaque latitude. +c angle zenithal calculé pour obtenir un +c Fit polynomial de l'ensoleillement moyen au sommet de l'atmosphere +c en moyenne annuelle. +c Spécifique de la terre. Utilisé pour les aqua planetes. +c +c Rque : Different de la routine angle en ce sens que zenang +c fournit des moyennes de pmu0 et non des valeurs +c instantanees, du coup frac prend toutes les valeurs +c entre 0 et 1. +c Date : premiere version le 13 decembre 1994 +c revu pour GCM le 30 septembre 1996 +c=============================================================== +c longi----INPUT : la longitude vraie de la terre dans son plan +c solaire a partir de l'equinoxe de printemps (degre) +c gmtime---INPUT : temps universel en fraction de jour +c pdtrad---INPUT : pas de temps du rayonnement (secondes) +c lat------INPUT : latitude en degres +c long-----INPUT : longitude en degres +c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad +c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad +c================================================================ +#include "YOMCST.h" +c================================================================ + logical cycle_diurne + real gmtime + real rlat(klon), rlon(klon), rmu0(klon), fract(klon) +c================================================================ + integer i + real gmtime1, gmtime2 + real pi_local + + + real rmu0m(klon),rmu0a(klon) +c + + pi_local = 4.0 * ATAN(1.0) + +c================================================================ +c Calcul de l'angle zenithal moyen sur la journee +c================================================================ + + DO i=1,klon + fract(i)=1. +! Calcule du flux moyen + IF (abs(rlat(i)).LE.28.75) THEN + rmu0m(i)=(210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365. + ELSEIF (abs(rlat(i)).LE.43.75) THEN + rmu0m(i)=(187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365. + ELSEIF (abs(rlat(i)).LE.71.25) THEN + rmu0m(i)=(162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365. + ELSE + rmu0m(i)=(172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365. + ENDIF + ENDDO + +c================================================================ +! Avec ou sans cycle diurne +c================================================================ + + IF (cycle_diurne) THEN + +! On redecompose flux au sommet suivant un cycle diurne idealise +! identique a toutes les latitudes. + + DO i=1,klon + rmu0a(i)=2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local) + rmu0(i)=rmu0a(i)*abs(sin(pi_local*gmtime+pi_local* + & rlon(i)/360.))-rmu0a(i)/sqrt(2.) + ENDDO + + DO i=1,klon + IF (rmu0(i).LE.0.) THEN + rmu0(i)=0. + fract(i)=0. + ELSE + fract(i)=1. + ENDIF + ENDDO + +! Affichage de l'angel zenitale +! print*,'************************************' +! print*,'************************************' +! print*,'************************************' +! print*,'latitude=',rlat(i),'longitude=',rlon(i) +! print*,'rmu0m=',rmu0m(i) +! print*,'rmu0a=',rmu0a(i) +! print*,'rmu0=',rmu0(i) + + ELSE + + DO i=1,klon + fract(i)=0.5 + rmu0(i)=rmu0m(i)*2. + ENDDO + + ENDIF + + RETURN + END + subroutine writelim + s (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice, + s phy_fter,phy_foce,phy_flic,phy_fsic) +c +!#include "dimensions.h" +!#include "dimphy.h" +#include "netcdf.inc" + + integer klon + REAL phy_nat(klon,360) + REAL phy_alb(klon,360) + REAL phy_sst(klon,360) + REAL phy_bil(klon,360) + REAL phy_rug(klon,360) + REAL phy_ice(klon,360) + REAL phy_fter(klon,360) + REAL phy_foce(klon,360) + REAL phy_flic(klon,360) + REAL phy_fsic(klon,360) + + INTEGER ierr + INTEGER dimfirst(3) + INTEGER dimlast(3) +c + INTEGER nid, ndim, ntim + INTEGER dims(2), debut(2), epais(2) + INTEGER id_tim + INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB + INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC + + PRINT*, 'Ecriture du fichier limit' +c + ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid) +c + ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30, + . "Fichier conditions aux limites") + ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim) + ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim) +c + dims(1) = ndim + dims(2) = ntim +c +ccc ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim) + ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim) + ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17, + . "Jour dans l annee") +ccc ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT) + ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT) + ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23, + . "Nature du sol (0,1,2,3)") +ccc ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST) + ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST) + ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35, + . "Temperature superficielle de la mer") +ccc ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS) + ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS) + ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32, + . "Reference flux de chaleur au sol") +ccc ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB) + ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB) + ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19, + . "Albedo a la surface") +ccc ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG) + ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG) + ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8, + . "Rugosite") + + ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER) + ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre") + ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE) + ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre") + ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC) + ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre") + ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC) + ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre") +c + ierr = NF_ENDDEF(nid) +c + DO k = 1, 360 +c + debut(1) = 1 + debut(2) = k + epais(1) = klon + epais(2) = 1 +c + print*,'Instant ',k +#ifdef NC_DOUBLE + print*,'NC DOUBLE' + ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais,phy_fter(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais,phy_foce(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais,phy_fsic(1,k)) + ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais,phy_flic(1,k)) +#else + print*,'NC PAS DOUBLE' + ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k)) + ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais,phy_fter(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais,phy_foce(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais,phy_fsic(1,k)) + ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais,phy_flic(1,k)) + +#endif +c + ENDDO +c + ierr = NF_CLOSE(nid) +c + return + end + + SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst) + use dimphy + IMPLICIT none +c + INTEGER nlon, type_profil, i, k, j + REAL :: rlatd(nlon), phy_sst(nlon, 360) + INTEGER imn, imx, amn, amx, kmn, kmx + INTEGER p, pplus, nlat_max + parameter (nlat_max=72) + REAL x_anom_sst(nlat_max) +c + if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua' + do i=1,360 +c phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2 + +c Rajout fbrlmd + + if(type_profil.EQ.1)then +c Méthode 1 "Control" faible plateau à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**2) +c PI/3=1.047197551 + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273. + endif + enddo + endif + if(type_profil.EQ.2)then +c Méthode 2 "Flat" fort plateau à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273. + endif + enddo + endif + + + if (type_profil.EQ.3) then +c Méthode 3 "Qobs" plateau réel à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273. + endif + enddo + endif + + if (type_profil.EQ.4) then +c Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273. + endif + enddo + endif + + if (type_profil.EQ.5) then +c Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+2.+0.5*27.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273.+2. + endif + + enddo + endif + + if(type_profil.EQ.6)then +c Méthode 6 "cst" valeur constante de SST + do j=1,klon + phy_sst(j,i)=288. + enddo + endif + + + if(type_profil.EQ.7)then +c Méthode 7 "cst" valeur constante de SST +2 + do j=1,klon + phy_sst(j,i)=288.+2. + enddo + endif + + p=0 + if(type_profil.EQ.8)then +c Méthode 8 profil anomalies SST du modèle couplé AR4 + do j=1,klon + if (rlatd(j).EQ.rlatd(j-1)) then + phy_sst(j,i)=273.+x_anom_sst(pplus) + & +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273.+x_anom_sst(pplus) + endif + else + p=p+1 + pplus=73-p + phy_sst(j,i)=273.+x_anom_sst(pplus) + & +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273.+x_anom_sst(pplus) + endif + write (*,*) rlatd(j),x_anom_sst(pplus),phy_sst(j,i) + endif + enddo + endif + + if (type_profil.EQ.9) then +c Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur + do j=1,klon + phy_sst(j,i)=273.-2.+0.5*27.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273.-2. + endif + enddo + endif + + + if (type_profil.EQ.10) then +c Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273.+4. + endif + enddo + endif + + if (type_profil.EQ.11) then +c Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273. + endif + enddo + endif + + if (type_profil.EQ.12) then +c Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur + do j=1,klon + phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273.+4. + endif + enddo + endif + + if (type_profil.EQ.13) then +c Méthode 13 "Qmax" plateau réel à l'Equateur augmenté ! + do j=1,klon + phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273. + endif + enddo + endif + + if (type_profil.EQ.14) then +c Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K ! + do j=1,klon + phy_sst(j,i)=273.+2.+0.5*29.*(2-sin(1.5*rlatd(j))**2 + & -sin(1.5*rlatd(j))**4) + if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then + phy_sst(j,i)=273. + endif + enddo + endif + + enddo + +cIM beg : verif profil SST: phy_sst + amn=MIN(phy_sst(1,1),1000.) + amx=MAX(phy_sst(1,1),-1000.) + DO k=1, 360 + DO i=2, nlon + IF(phy_sst(i,k).LT.amn) THEN + amn=phy_sst(i,k) + imn=i + kmn=k + ENDIF + IF(phy_sst(i,k).GT.amx) THEN + amx=phy_sst(i,k) + imx=i + kmx=k + ENDIF + ENDDO + ENDDO +c + PRINT*,' debut avant writelim min max phy_sst',imn,kmn,amn, + & imx,kmx,amx +cIM end : verif profil SST: phy_sst + + return + end Index: LMDZ5/trunk/libf/phylmd/physiq.F =================================================================== --- LMDZ5/trunk/libf/phylmd/physiq.F (revision 1527) +++ LMDZ5/trunk/libf/phylmd/physiq.F (revision 1529) @@ -1881,10 +1881,22 @@ & write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Calcul de l'ensoleillement : +! ============================ +! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur +! l'annee a partir d'une formule analytique. +! Cet ensoleillement est symmétrique autour de l'équateur et +! non nul aux poles. + IF (abs(solarlong0-1000.)<1.e-4) then + call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract) + ELSE ! Avec ou sans cycle diurne - IF (cycle_diurne) THEN - zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s) - CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract) - ELSE - CALL angle(zlongi, rlat, fract, rmu0) + IF (cycle_diurne) THEN + zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s) + CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract) + ELSE + CALL angle(zlongi, rlat, fract, rmu0) + ENDIF ENDIF