! Routines complementaires pour la physique planetaire. MODULE phyaqua_mod CONTAINS 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, only : rlatd,rlond USE dimphy, only : klon 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, only : dayref,nday,iphysiq USE indice_sol_mod USE IOIPSL IMPLICIT NONE #include "dimensions.h" ! #include "dimphy.h" ! #include "YOMCST.h" #include "comconst.h" #include "clesphys.h" #include "dimsoil.h" #include "temps.h" integer,intent(in) :: nlon,iflag_phys !IM ajout latfi, lonfi real,intent(in) :: lonfi(nlon),latfi(nlon) 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 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=0 ! read ozone climatology ! intermediate variables to use getin integer :: nbapp_rad_omp real :: co2_ppm_omp,solaire_omp logical :: alb_ocean_omp real :: rugos_omp !------------------------------------------------------------------------- ! 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*,'iniaqua:type_aqua, type_profil',type_aqua, type_profil if (klon.ne.nlon) then write(*,*)"iniaqua: klon=",klon," nlon=",nlon stop'probleme de dimensions dans iniaqua' endif call phys_state_var_init(read_climoz) read_climoz=0 day0=217. day=day0 it=0 time=0. !IM ajout latfi, lonfi rlatd=latfi rlond=lonfi rlat=rlatd*180./pi rlon=rlond*180./pi !----------------------------------------------------------------------- ! initialisations de la physique !----------------------------------------------------------------------- day_ini=dayref day_end=day_ini+nday airefi=1. zcufi=1. zcvfi=1. !$OMP MASTER nbapp_rad_omp=24 CALL getin('nbapp_rad',nbapp_rad_omp) !$OMP END MASTER !$OMP BARRIER nbapp_rad=nbapp_rad_omp !--------------------------------------------------------------------- !c Creation des conditions aux limites: !c ------------------------------------ ! Initialisations des constantes ! Ajouter les manquants dans planete.def... (albedo etc) !$OMP MASTER co2_ppm_omp=348. CALL getin('co2_ppm',co2_ppm_omp) solaire_omp=1365. CALL getin('solaire',solaire_omp) ! CALL getin('albedo',albedo) ! albedo is set below, depending on type_aqua alb_ocean_omp=.true. CALL getin('alb_ocean',alb_ocean_omp) !$OMP END MASTER !$OMP BARRIER co2_ppm=co2_ppm_omp solaire=solaire_omp alb_ocean=alb_ocean_omp radsol=0. qsol_f=10. !c Conditions aux limites: !c ----------------------- qsol(:) = qsol_f rugsrel = 0.0 ! (rugsrel = rugoro) rugoro = 0.0 u_ancien = 0.0 v_ancien = 0.0 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 !$OMP MASTER rugos_omp=rugos CALL getin('rugos',rugos_omp) !$OMP END MASTER !$OMP BARRIER rugos=rugos_omp 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 !IM Terraplanete phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2 !IM 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 !IM calcul profil sst call profil_sst(nlon, rlatd, type_profil, phy_sst) call writelim & & (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice, & & 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*,'iniaqua: before phyredem' falb1=albedo falb2=albedo zmax0=0. f0=0. ema_work1=0. ema_work2=0. wake_deltat=0. wake_deltaq=0. wake_s=0. wake_cstar=0. wake_pe=0. wake_fip=0. fm_therm=0. entr_therm=0. detr_therm=0. CALL phyredem ("startphy.nc") print*,'iniaqua: after phyredem' call phys_state_var_end return end SUBROUTINE iniaqua !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 zenang_an !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine writelim & & (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice, & & phy_fter,phy_foce,phy_flic,phy_fsic) !c use mod_phys_lmdz_para, only: is_mpi_root,is_omp_root use mod_grid_phy_lmdz, only : klon_glo use mod_phys_lmdz_transfert_para, only : gather !L. Fita July 2013 use netcdf !#include "dimensions.h" !#include "dimphy.h" integer,intent(in) :: klon real,intent(in) :: phy_nat(klon,360) real,intent(in) :: phy_alb(klon,360) real,intent(in) :: phy_sst(klon,360) real,intent(in) :: phy_bil(klon,360) real,intent(in) :: phy_rug(klon,360) real,intent(in) :: phy_ice(klon,360) real,intent(in) :: phy_fter(klon,360) real,intent(in) :: phy_foce(klon,360) real,intent(in) :: phy_flic(klon,360) real,intent(in) :: phy_fsic(klon,360) real :: phy_glo(klon_glo,360) ! temporary variable, to store phy_***(:) ! on the whole physics grid 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 if (is_mpi_root.and.is_omp_root) then PRINT*, 'writelim: 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, "points_physiques", klon_glo, 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 ! write the 'times' do k=1,360 #ifdef NC_DOUBLE ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k)) #else ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k)) #endif enddo endif ! of if (is_mpi_root.and.is_omp_root) ! write the fields, after having collected them on master call gather(phy_nat,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_NAT,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_NAT,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_nat" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_sst,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_SST,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_SST,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_sst" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_bil,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_BILS,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_BILS,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_bil" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_alb,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_ALB,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_ALB,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_alb" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_rug,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_RUG,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_RUG,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_rug" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_fter,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_FTER,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_FTER,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_fter" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_foce,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_FOCE,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_FOCE,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_foce" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_fsic,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_FSIC,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_FSIC,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_fsic" write(*,*) NF_STRERROR(ierr) endif endif call gather(phy_flic,phy_glo) if (is_mpi_root.and.is_omp_root) then #ifdef NC_DOUBLE ierr=NF_PUT_VAR_DOUBLE(nid,id_FLIC,phy_glo) #else ierr=NF_PUT_VAR_REAL(nid,id_FLIC,phy_glo) #endif if(ierr.ne.NF_NOERR) then write(*,*) "writelim error with phy_flic" write(*,*) NF_STRERROR(ierr) endif endif ! close file: if (is_mpi_root.and.is_omp_root) then ierr = NF_CLOSE(nid) endif end SUBROUTINE writelim !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 !IM 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 !IM end : verif profil SST: phy_sst return end SUBROUTINE profil_sst END MODULE phyaqua_mod