MODULE phyaqua_mod ! Routines complementaires pour la physique planetaire. IMPLICIT NONE CONTAINS SUBROUTINE iniaqua(nlon, 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 dimphy, ONLY: klon USE geometry_mod, ONLY : latitude 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 time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, & day_ini,day_end USE indice_sol_mod USE nrtype, ONLY: pi USE ioipsl #ifdef ISO USE infotrac_phy, ONLY: niso #endif IMPLICIT NONE include "YOMCST.h" include "clesphys.h" include "dimsoil.h" INTEGER, INTENT (IN) :: nlon, iflag_phys ! IM ajout latfi, lonfi ! REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon) INTEGER type_profil, type_aqua ! Ajouts initialisation des surfaces REAL :: run_off_lic_0(nlon) REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf) REAL :: tsoil(nlon, nsoilmx, nbsrf) REAL :: tslab(nlon), seaice(nlon) REAL fder(nlon) #ifdef ISO REAL :: xtrun_off_lic_0(niso,nlon) REAL :: xtsolsrf(niso,nlon, nbsrf), xtsnsrf(niso,nlon, nbsrf) REAL :: Rland_ice(niso,nlon) #endif ! Arguments : ! ----------- ! 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 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 (need to be "save" to be shared by ! all threads) INTEGER, SAVE :: nbapp_rad_omp REAL, SAVE :: co2_ppm_omp, solaire_omp LOGICAL, SAVE :: alb_ocean_omp REAL, SAVE :: 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) ! ----------------------------------------------------------------------- ! dynamial tendencies : ! --------------------- 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/=nlon) THEN WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon STOP 'probleme de dimensions dans iniaqua' END IF 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 = day_ref day_end = day_ini + ndays ! 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 ! --------------------------------------------------------------------- ! Creation des conditions aux limites: ! ------------------------------------ ! 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 WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm solaire = solaire_omp WRITE (*, *) 'iniaqua: solaire=', solaire alb_ocean = alb_ocean_omp WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean radsol = 0. qsol_f = 10. ! Conditions aux limites: ! ----------------------- 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. END IF !$OMP MASTER rugos_omp = rugos CALL getin('rugos', rugos_omp) !$OMP END MASTER !$OMP BARRIER rugos = rugos_omp WRITE (*, *) 'iniaqua: rugos=', rugos zmasq(:) = pctsrf(:, is_ter) ! 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 (?) END DO ! IM calcul profil sst CALL profil_sst(nlon, latitude, 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) ! --------------------------------------------------------------------- ! Ecriture de l'etat initial: ! --------------------------- ! Ecriture etat initial physique timestep = pdtphys radpas = nint(rday/timestep/float(nbapp_rad)) DO i = 1, longcles clesphy0(i) = 0. END DO clesphy0(1) = float(iflag_con) clesphy0(2) = float(nbapp_rad) ! 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 ! ======================================================================= ! Profils initiaux ! ======================================================================= ! 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) END DO 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. fevap = 0. #ifdef ISO xtrun_off_lic_0 = 0. #endif ! Initialisations necessaires avant phyredem type_ocean = 'force' CALL fonte_neige_init(run_off_lic_0 & #ifdef ISO & ,xtrun_off_lic_0 & #endif & ) qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface snsrf(:, :) = 0. ! couverture de neige des sous surface z0m(:, :) = rugos ! couverture de neige des sous surface z0h=z0m #ifdef ISO xtsolsrf(:, :,:) = qsol(1) xtsnsrf(:, :,:) = 0. #endif CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil & #ifdef ISO & ,xtsnsrf,Rland_ice & #endif & ) PRINT *, 'iniaqua: before phyredem' falb1 = albedo falb2 = albedo zmax0 = 0. f0 = 0. sig1 = 0. w01 = 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 ! ==================================================================== ! ==================================================================== SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract) USE dimphy IMPLICIT NONE ! ==================================================================== ! ============================================================= ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract) ! Auteur : A. Campoy et F. Hourdin ! Objet : calculer les valeurs moyennes du cos de l'angle zenithal ! et l'ensoleillement moyen entre gmtime1 et gmtime2 ! connaissant la declinaison, la latitude et la longitude. ! Dans cette version particuliere, on calcule le rayonnement ! moyen sur l'année à chaque latitude. ! angle zenithal calculé pour obtenir un ! Fit polynomial de l'ensoleillement moyen au sommet de l'atmosphere ! en moyenne annuelle. ! Spécifique de la terre. Utilisé pour les aqua planetes. ! Rque : Different de la routine angle en ce sens que zenang ! fournit des moyennes de pmu0 et non des valeurs ! instantanees, du coup frac prend toutes les valeurs ! entre 0 et 1. ! Date : premiere version le 13 decembre 1994 ! revu pour GCM le 30 septembre 1996 ! =============================================================== ! longi----INPUT : la longitude vraie de la terre dans son plan ! solaire a partir de l'equinoxe de printemps (degre) ! gmtime---INPUT : temps universel en fraction de jour ! pdtrad---INPUT : pas de temps du rayonnement (secondes) ! lat------INPUT : latitude en degres ! long-----INPUT : longitude en degres ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad ! ================================================================ include "YOMCST.h" ! ================================================================ LOGICAL cycle_diurne REAL gmtime REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon) ! ================================================================ INTEGER i REAL gmtime1, gmtime2 REAL pi_local REAL rmu0m(klon), rmu0a(klon) pi_local = 4.0*atan(1.0) ! ================================================================ ! Calcul de l'angle zenithal moyen sur la journee ! ================================================================ DO i = 1, klon fract(i) = 1. ! Calcule du flux moyen IF (abs(rlat(i))<=28.75) THEN rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365. ELSE IF (abs(rlat(i))<=43.75) THEN rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365. ELSE IF (abs(rlat(i))<=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. END IF END DO ! ================================================================ ! Avec ou sans cycle diurne ! ================================================================ 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.) END DO DO i = 1, klon IF (rmu0(i)<=0.) THEN rmu0(i) = 0. fract(i) = 0. ELSE fract(i) = 1. END IF END DO ! 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. END DO END IF 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) 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 IMPLICIT NONE include "netcdf.inc" 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 :: k INTEGER ierr INTEGER dimfirst(3) INTEGER dimlast(3) 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' ierr = nf_create('limit.nc', nf_clobber, nid) 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) dims(1) = ndim dims(2) = ntim #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim) #else ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim) #endif ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat) #else ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat) #endif ierr = nf_put_att_text(nid, id_nat, 'title', 23, & 'Nature du sol (0,1,2,3)') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst) #else ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst) #endif ierr = nf_put_att_text(nid, id_sst, 'title', 35, & 'Temperature superficielle de la mer') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils) #else ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils) #endif ierr = nf_put_att_text(nid, id_bils, 'title', 32, & 'Reference flux de chaleur au sol') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb) #else ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb) #endif ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug) #else ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug) #endif ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter) #else ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter) #endif ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce) #else ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce) #endif ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic) #else ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic) #endif ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice') #ifdef NC_DOUBLE ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic) #else ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic) #endif ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice') ierr = nf_enddef(nid) IF (ierr/=nf_noerr) THEN WRITE (*, *) 'writelim error: failed to end define mode' WRITE (*, *) nf_strerror(ierr) END IF ! 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 IF (ierr/=nf_noerr) THEN WRITE (*, *) 'writelim error with temps(k),k=', k WRITE (*, *) nf_strerror(ierr) END IF END DO END IF ! 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_nat' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_sst' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_bil' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_alb' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_rug' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_fter' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_foce' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_fsic' WRITE (*, *) nf_strerror(ierr) END IF END IF 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/=nf_noerr) THEN WRITE (*, *) 'writelim error with phy_flic' WRITE (*, *) nf_strerror(ierr) END IF END IF ! close file: IF (is_mpi_root .AND. is_omp_root) THEN ierr = nf_close(nid) END IF END SUBROUTINE writelim ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst) USE dimphy IMPLICIT NONE 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) IF (klon/=nlon) STOP 'probleme de dimensions dans iniaqua' WRITE (*, *) ' profil_sst: type_profil=', type_profil DO i = 1, 360 ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2 ! Rajout fbrlmd IF (type_profil==1) THEN ! 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) ! PI/3=1.047197551 IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. END IF END DO END IF IF (type_profil==2) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. END IF END DO END IF IF (type_profil==3) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. END IF END DO END IF IF (type_profil==4) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. END IF END DO END IF IF (type_profil==5) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. + 2. END IF END DO END IF IF (type_profil==6) THEN ! Méthode 6 "cst" valeur constante de SST DO j = 1, klon phy_sst(j, i) = 288. END DO END IF IF (type_profil==7) THEN ! Méthode 7 "cst" valeur constante de SST +2 DO j = 1, klon phy_sst(j, i) = 288. + 2. END DO END IF p = 0 IF (type_profil==8) THEN ! Méthode 8 profil anomalies SST du modèle couplé AR4 DO j = 1, klon IF (rlatd(j)==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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. + x_anom_sst(pplus) END IF 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. + x_anom_sst(pplus) END IF WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i) END IF END DO END IF IF (type_profil==9) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. - 2. END IF END DO END IF IF (type_profil==10) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. + 4. END IF END DO END IF IF (type_profil==11) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. END IF END DO END IF IF (type_profil==12) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. + 4. END IF END DO END IF IF (type_profil==13) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. END IF END DO END IF IF (type_profil==14) THEN ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN phy_sst(j, i) = 273. END IF END DO END IF if (type_profil.EQ.20) then print*,'Profile SST 20' ! Méthode 13 "Qmax2K" plateau réel �| l'Equateur augmenté +2K do j=1,klon phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2) enddo endif if (type_profil.EQ.21) then print*,'Profile SST 21' ! Méthode 13 "Qmax2K" plateau réel �| l'Equateur augmenté +2K do j=1,klon phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2) enddo endif END DO ! IM beg : verif profil SST: phy_sst amn = min(phy_sst(1,1), 1000.) amx = max(phy_sst(1,1), -1000.) imn = 1 kmn = 1 imx = 1 kmx = 1 DO k = 1, 360 DO i = 2, nlon IF (phy_sst(i,k)amx) THEN amx = phy_sst(i, k) imx = i kmx = k END IF END DO END DO PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx ! IM end : verif profil SST: phy_sst RETURN END SUBROUTINE profil_sst END MODULE phyaqua_mod