Context Navigation

phyaqua_mod.F90 @ 5119

Last change on this file since 5119 was 5117, checked in by abarral, 4 months ago
rename modules properly lmdz_* move some unused files to obsolete/ (lint) uppercase fortran keywords
File size: 30.5 KB

Rev	Line
[5099]	1
[3927]	2	! $Id: phyaqua_mod.F90 3579 2019-10-09 13:11:07Z fairhead $
[5099]	3
[3927]	4	MODULE phyaqua_mod
	5	! Routines complementaires pour la physique planetaire.
	6	IMPLICIT NONE
	7
	8	CONTAINS
	9
	10	SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
	11
	12	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	13	! Creation d'un etat initial et de conditions aux limites
	14	! (resp startphy.nc et limit.nc) pour des configurations idealisees
	15	! du modele LMDZ dans sa version terrestre.
	16	! iflag_phys est un parametre qui controle
	17	! iflag_phys = N
	18	! de 100 a 199 : aqua planetes avec SST forcees
	19	! N-100 determine le type de SSTs
	20	! de 200 a 299 : terra planetes avec Ts calcule
	21
	22	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	23
	24	USE dimphy, ONLY: klon
[5112]	25	USE lmdz_geometry, ONLY: latitude
[3927]	26	USE surface_data, ONLY: type_ocean, ok_veget
	27	USE pbl_surface_mod, ONLY: pbl_surface_init
	28	USE fonte_neige_mod, ONLY: fonte_neige_init
	29	#ifdef ISO
[5101]	30	USE fonte_neige_mod, ONLY: fonte_neige_init_iso
	31	USE pbl_surface_mod, ONLY: pbl_surface_init_iso
[3927]	32	#endif
	33	USE phys_state_var_mod
	34	USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
	35	day_ini,day_end
	36	USE indice_sol_mod
[5117]	37	USE lmdz_physical_constants, ONLY: pi
[3927]	38	! USE ioipsl
[5110]	39	USE lmdz_phys_para, ONLY: is_master
	40	USE lmdz_phys_transfert_para, ONLY: bcast
	41	USE lmdz_grid_phy
[5112]	42	USE lmdz_ioipsl_getin_p, ONLY: getin_p
[3927]	43	USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
	44	#ifdef ISO
	45	USE infotrac_phy, ONLY: niso
	46	#endif
	47
	48	IMPLICIT NONE
	49
	50	include "YOMCST.h"
	51	include "clesphys.h"
	52	include "dimsoil.h"
	53
	54	INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
	55	! IM ajout latfi, lonfi
	56	! REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
	57
	58	INTEGER type_profil, type_aqua
	59
	60	! Ajouts initialisation des surfaces
	61	REAL :: run_off_lic_0(nlon)
	62	REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
	63	REAL :: tsoil(nlon, nsoilmx, nbsrf)
	64	REAL :: tslab(nlon), seaice(nlon)
	65	REAL fder(nlon)
	66
	67	#ifdef ISO
	68	REAL :: xtrun_off_lic_0(niso,nlon)
	69	REAL :: xtsolsrf(niso,nlon, nbsrf), xtsnsrf(niso,nlon, nbsrf)
	70	REAL :: Rland_ice(niso,nlon)
	71	#endif
	72
	73
	74
	75	! Arguments :
	76	! -----------
	77
	78	! integer radpas
	79	INTEGER it, unit, i, k, itap
	80
	81	REAL rugos, albedo
	82	REAL tsurf
	83	REAL time, timestep, day, day0
	84	REAL qsol_f
	85	REAL rugsrel(nlon)
	86	LOGICAL alb_ocean
	87
	88	CHARACTER *80 ans, file_forctl, file_fordat, file_start
	89	CHARACTER *100 file, var
	90	CHARACTER *2 cnbl
	91
	92	REAL phy_nat(nlon, year_len)
	93	REAL phy_alb(nlon, year_len)
	94	REAL phy_sst(nlon, year_len)
	95	REAL phy_bil(nlon, year_len)
	96	REAL phy_rug(nlon, year_len)
	97	REAL phy_ice(nlon, year_len)
	98	REAL phy_fter(nlon, year_len)
	99	REAL phy_foce(nlon, year_len)
	100	REAL phy_fsic(nlon, year_len)
	101	REAL phy_flic(nlon, year_len)
	102
	103	INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
	104	!$OMP THREADPRIVATE(read_climoz)
	105
	106	! -------------------------------------------------------------------------
	107	! declaration pour l'appel a phyredem
	108	! -------------------------------------------------------------------------
	109
	110	! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
	111	REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
	112	! real pbl_tke(nlon,llm,nbsrf)
	113	! real rain_fall(nlon),snow_fall(nlon)
	114	! real solsw(nlon), sollw(nlon),radsol(nlon)
	115	! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
	116	! real ratqs(nlon,llm)
	117	! real clwcon(nlon,llm)
	118
	119	INTEGER longcles
	120	PARAMETER (longcles=20)
	121	REAL clesphy0(longcles)
	122
	123
	124	! -----------------------------------------------------------------------
	125	! dynamial tendencies :
	126	! ---------------------
	127
	128	INTEGER l, ierr, aslun
	129
	130	REAL paire
	131
	132	! Local
	133	CHARACTER (LEN=20) :: modname='phyaqua'
	134	CHARACTER (LEN=80) :: abort_message
	135
	136
	137	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	138	! INITIALISATIONS
	139	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	140
	141	! -----------------------------------------------------------------------
	142	! Initialisations des constantes
	143	! -------------------------------
	144
	145	!IF (calend .EQ. "earth_360d") Then
	146	year_len_phy = year_len
	147	!END IF
	148
[5117]	149	IF (year_len/=360) THEN
[3927]	150	write (,) year_len
	151	write (,) 'iniaqua: 360 day calendar is required !'
	152	stop
	153	endif
	154
	155	type_aqua = iflag_phys/100
	156	type_profil = iflag_phys - type_aqua*100
	157	PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
	158
	159	IF (klon/=nlon) THEN
	160	WRITE (, ) 'iniaqua: klon=', klon, ' nlon=', nlon
	161	abort_message= 'probleme de dimensions dans iniaqua'
	162	CALL abort_physic(modname,abort_message,1)
	163	END IF
	164	CALL phys_state_var_init(read_climoz)
	165
	166
	167	read_climoz = 0
	168	day0 = 217.
	169	day = day0
	170	it = 0
	171	time = 0.
	172
	173	! -----------------------------------------------------------------------
	174	! initialisations de la physique
	175	! -----------------------------------------------------------------------
	176
	177	day_ini = day_ref
	178	day_end = day_ini + ndays
	179
	180	nbapp_rad = 24
	181	CALL getin_p('nbapp_rad', nbapp_rad)
	182
	183	! ---------------------------------------------------------------------
	184	! Creation des conditions aux limites:
	185	! ------------------------------------
	186	! Initialisations des constantes
	187	! Ajouter les manquants dans planete.def... (albedo etc)
	188	co2_ppm = 348.
	189	CALL getin_p('co2_ppm', co2_ppm)
	190
	191	solaire = 1365.
	192	CALL getin_p('solaire', solaire)
	193
	194	! CALL getin('albedo',albedo) ! albedo is set below, depending on
	195	! type_aqua
	196	alb_ocean = .TRUE.
	197	CALL getin_p('alb_ocean', alb_ocean)
	198
	199	WRITE (, ) 'iniaqua: co2_ppm=', co2_ppm
	200	WRITE (, ) 'iniaqua: solaire=', solaire
	201	WRITE (, ) 'iniaqua: alb_ocean=', alb_ocean
	202
	203	radsol = 0.
	204	qsol_f = 10.
	205
	206	! Conditions aux limites:
	207	! -----------------------
	208
	209	qsol(:) = qsol_f
	210	rugsrel = 0.0 ! (rugsrel = rugoro)
	211	rugoro = 0.0
	212	u_ancien = 0.0
	213	v_ancien = 0.0
	214	agesno = 50.0
	215	! Relief plat
	216	zmea = 0.
	217	zstd = 0.
	218	zsig = 0.
	219	zgam = 0.
	220	zthe = 0.
	221	zpic = 0.
	222	zval = 0.
	223
	224	! Une seule surface
	225	pctsrf = 0.
	226	IF (type_aqua==1) THEN
	227	rugos = 1.E-4
	228	albedo = 0.19
	229	pctsrf(:, is_oce) = 1.
	230	ELSE IF (type_aqua==2) THEN
	231	rugos = 0.03
	232	albedo = 0.1
	233	pctsrf(:, is_ter) = 1.
	234	END IF
	235
	236	CALL getin_p('rugos', rugos)
	237
	238	WRITE (, ) 'iniaqua: rugos=', rugos
	239	zmasq(:) = pctsrf(:, is_ter)
	240
	241	! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
	242	! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
	243
	244	! Si alb_ocean on calcule un albedo oceanique moyen
[5116]	245	! if (alb_ocean) THEN
[3927]	246	! Voir pourquoi on avait ca.
	247	! CALL ini_alb_oce(phy_alb)
	248	! else
	249	phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
[5117]	250	! END IF !alb_ocean
[3927]	251
	252	DO i = 1, year_len
	253	! IM Terraplanete phy_sst(:,i) = 260.+50.cos(rlatd(:))*2
	254	! IM ajout calcul profil sst selon le cas considere (cf. FBr)
	255
	256	phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
	257	phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
	258	phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
	259	phy_ice(:, i) = 0.0 ! fraction de glace (?)
	260	phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
	261	phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
	262	phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
	263	phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
	264	END DO
	265	! IM calcul profil sst
	266	CALL profil_sst(nlon, latitude, type_profil, phy_sst)
	267
	268	IF (grid_type==unstructured) THEN
	269	CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
	270	phy_fter, phy_foce, phy_flic, phy_fsic)
	271	ELSE
	272
	273	CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
	274	phy_fter, phy_foce, phy_flic, phy_fsic)
	275	ENDIF
	276
	277	! ---------------------------------------------------------------------
	278	! Ecriture de l'etat initial:
	279	! ---------------------------
	280
	281
	282	! Ecriture etat initial physique
	283
	284	timestep = pdtphys
	285	radpas = nint(rday/timestep/float(nbapp_rad))
	286
	287	DO i = 1, longcles
	288	clesphy0(i) = 0.
	289	END DO
	290	clesphy0(1) = float(iflag_con)
	291	clesphy0(2) = float(nbapp_rad)
	292	! IF( cycle_diurne ) clesphy0(3) = 1.
	293	clesphy0(3) = 1. ! cycle_diurne
	294	clesphy0(4) = 1. ! soil_model
	295	clesphy0(5) = 1. ! new_oliq
	296	clesphy0(6) = 0. ! ok_orodr
	297	clesphy0(7) = 0. ! ok_orolf
	298	clesphy0(8) = 0. ! ok_limitvrai
	299
	300
	301	! =======================================================================
	302	! Profils initiaux
	303	! =======================================================================
	304
	305	! On initialise les temperatures de surfaces comme les sst
	306	DO i = 1, nlon
	307	ftsol(i, :) = phy_sst(i, 1)
	308	tsoil(i, :, :) = phy_sst(i, 1)
	309	tslab(i) = phy_sst(i, 1)
	310	END DO
	311
	312	falbe(:, :) = albedo
	313	falblw(:, :) = albedo
	314	rain_fall(:) = 0.
	315	snow_fall(:) = 0.
	316	solsw(:) = 0.
	317	sollw(:) = 0.
	318	radsol(:) = 0.
	319
	320	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	321	! intialisation bidon mais pas grave
	322	t_ancien(:, :) = 0.
	323	q_ancien(:, :) = 0.
	324	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	325	rnebcon = 0.
	326	ratqs = 0.
	327	clwcon = 0.
	328	pbl_tke = 1.E-8
	329
	330	! variables supplementaires pour appel a plb_surface_init
	331	fder(:) = 0.
	332	seaice(:) = 0.
	333	run_off_lic_0 = 0.
	334	fevap = 0.
	335	#ifdef ISO
	336	xtrun_off_lic_0 = 0.
	337	#endif
	338
	339
	340	! Initialisations necessaires avant phyredem
	341	type_ocean = 'force'
	342	CALL fonte_neige_init(run_off_lic_0)
	343	#ifdef ISO
	344	CALL fonte_neige_init_iso(xtrun_off_lic_0)
	345	#endif
	346	qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
	347	snsrf(:, :) = 0. ! couverture de neige des sous surface
	348	z0m(:, :) = rugos ! couverture de neige des sous surface
	349	z0h=z0m
	350
	351
	352	CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
	353	#ifdef ISO
	354	CALL pbl_surface_init_iso(xtsnsrf,Rland_ice)
	355	#endif
	356
	357	PRINT *, 'iniaqua: before phyredem'
	358
	359	pbl_tke(:,:,:) = 1.e-8
	360	falb1 = albedo
	361	falb2 = albedo
	362	zmax0 = 0.
	363	f0 = 0.
	364	sig1 = 0.
	365	w01 = 0.
	366	wake_deltat = 0.
	367	wake_deltaq = 0.
	368	wake_s = 0.
	369	wake_dens = 0.
	370	wake_cstar = 0.
	371	wake_pe = 0.
	372	wake_fip = 0.
	373	fm_therm = 0.
	374	entr_therm = 0.
	375	detr_therm = 0.
	376	ale_bl = 0.
	377	ale_bl_trig =0.
	378	alp_bl =0.
	379	treedrg(:,:,:)=0.
	380
	381	u10m = 0.
	382	v10m = 0.
	383
	384	ql_ancien = 0.
	385	qs_ancien = 0.
	386	u_ancien = 0.
	387	v_ancien = 0.
	388	prw_ancien = 0.
	389	prlw_ancien = 0.
	390	prsw_ancien = 0.
	391
	392	ale_wake = 0.
	393	ale_bl_stat = 0.
	394
	395
	396	!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
	397	! falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
	398	! falb_dir(:,is_oce,:)=0.5; falb_dir(:,is_sic,:)=0.6
	399	falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
	400	falb_dir(:,:,is_oce)=0.5; falb_dir(:,:,is_sic)=0.6
	401
	402	!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
	403	!ym probably the uninitialized value was 0 for standard (regular grid) case
	404	falb_dif(:,:,:)=0
	405
	406
	407	CALL phyredem('startphy.nc')
	408
	409	PRINT *, 'iniaqua: after phyredem'
	410	CALL phys_state_var_end
	411
[5105]	412
[3927]	413	END SUBROUTINE iniaqua
	414
	415
	416	! ====================================================================
	417	! ====================================================================
	418	SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
	419	USE dimphy
	420	IMPLICIT NONE
	421	! ====================================================================
	422	! =============================================================
	423	! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
	424	! Auteur : A. Campoy et F. Hourdin
	425	! Objet : calculer les valeurs moyennes du cos de l'angle zenithal
	426	! et l'ensoleillement moyen entre gmtime1 et gmtime2
	427	! connaissant la declinaison, la latitude et la longitude.
	428
	429	! Dans cette version particuliere, on calcule le rayonnement
	430	! moyen sur l'année à chaque latitude.
	431	! angle zenithal calculé pour obtenir un
	432	! Fit polynomial de l'ensoleillement moyen au sommet de l'atmosphere
	433	! en moyenne annuelle.
	434	! Spécifique de la terre. Utilisé pour les aqua planetes.
	435
	436	! Rque : Different de la routine angle en ce sens que zenang
	437	! fournit des moyennes de pmu0 et non des valeurs
	438	! instantanees, du coup frac prend toutes les valeurs
	439	! entre 0 et 1.
	440	! Date : premiere version le 13 decembre 1994
	441	! revu pour GCM le 30 septembre 1996
	442	! ===============================================================
	443	! longi----INPUT : la longitude vraie de la terre dans son plan
	444	! solaire a partir de l'equinoxe de printemps (degre)
	445	! gmtime---INPUT : temps universel en fraction de jour
	446	! pdtrad---INPUT : pas de temps du rayonnement (secondes)
	447	! lat------INPUT : latitude en degres
	448	! long-----INPUT : longitude en degres
	449	! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
	450	! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
	451	! ================================================================
	452	include "YOMCST.h"
	453	! ================================================================
	454	LOGICAL cycle_diurne
	455	REAL gmtime
	456	REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
	457	! ================================================================
	458	INTEGER i
	459	REAL gmtime1, gmtime2
	460	REAL pi_local
	461
	462
	463	REAL rmu0m(klon), rmu0a(klon)
	464
	465
	466	pi_local = 4.0*atan(1.0)
	467
	468	! ================================================================
	469	! Calcul de l'angle zenithal moyen sur la journee
	470	! ================================================================
	471
	472	DO i = 1, klon
	473	fract(i) = 1.
	474	! Calcule du flux moyen
	475	IF (abs(rlat(i))<=28.75) THEN
	476	rmu0m(i) = (210.1924+206.6059cos(0.0174533rlat(i))**2)/1365.
	477	ELSE IF (abs(rlat(i))<=43.75) THEN
	478	rmu0m(i) = (187.4562+236.1853cos(0.0174533rlat(i))**2)/1365.
	479	ELSE IF (abs(rlat(i))<=71.25) THEN
	480	rmu0m(i) = (162.4439+284.1192cos(0.0174533rlat(i))**2)/1365.
	481	ELSE
	482	rmu0m(i) = (172.8125+183.7673cos(0.0174533rlat(i))**2)/1365.
	483	END IF
	484	END DO
	485
	486	! ================================================================
	487	! Avec ou sans cycle diurne
	488	! ================================================================
	489
	490	IF (cycle_diurne) THEN
	491
	492	! On redecompose flux au sommet suivant un cycle diurne idealise
	493	! identique a toutes les latitudes.
	494
	495	DO i = 1, klon
	496	rmu0a(i) = 2.rmu0m(i)sqrt(2.)*pi_local/(4.-pi_local)
	497	rmu0(i) = rmu0a(i)abs(sin(pi_localgmtime+pi_local*rlon(i)/360.)) - &
	498	rmu0a(i)/sqrt(2.)
	499	END DO
	500
	501	DO i = 1, klon
	502	IF (rmu0(i)<=0.) THEN
	503	rmu0(i) = 0.
	504	fract(i) = 0.
	505	ELSE
	506	fract(i) = 1.
	507	END IF
	508	END DO
	509
	510	! Affichage de l'angel zenitale
[5103]	511	! PRINT,'***********************************'
	512	! PRINT,'***********************************'
	513	! PRINT,'***********************************'
	514	! PRINT*,'latitude=',rlat(i),'longitude=',rlon(i)
	515	! PRINT*,'rmu0m=',rmu0m(i)
	516	! PRINT*,'rmu0a=',rmu0a(i)
	517	! PRINT*,'rmu0=',rmu0(i)
[3927]	518
	519	ELSE
	520
	521	DO i = 1, klon
	522	fract(i) = 0.5
	523	rmu0(i) = rmu0m(i)*2.
	524	END DO
	525
	526	END IF
	527
[5105]	528
[3927]	529	END SUBROUTINE zenang_an
	530
	531	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	532
	533	SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
	534	phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
	535
[5110]	536	USE lmdz_phys_para, ONLY: is_omp_master, klon_mpi
	537	USE lmdz_phys_transfert_para, ONLY: gather_omp
[4619]	538	USE lmdz_xios
[3927]	539	IMPLICIT NONE
	540
	541	INTEGER, INTENT (IN) :: klon
	542	REAL, INTENT (IN) :: phy_nat(klon, 360)
	543	REAL, INTENT (IN) :: phy_alb(klon, 360)
	544	REAL, INTENT (IN) :: phy_sst(klon, 360)
	545	REAL, INTENT (IN) :: phy_bil(klon, 360)
	546	REAL, INTENT (IN) :: phy_rug(klon, 360)
	547	REAL, INTENT (IN) :: phy_ice(klon, 360)
	548	REAL, INTENT (IN) :: phy_fter(klon, 360)
	549	REAL, INTENT (IN) :: phy_foce(klon, 360)
	550	REAL, INTENT (IN) :: phy_flic(klon, 360)
	551	REAL, INTENT (IN) :: phy_fsic(klon, 360)
	552
	553	REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
	554	! on the whole physics grid
	555
	556	PRINT *, 'writelim: Ecriture du fichier limit'
	557
	558	CALL gather_omp(phy_foce, phy_mpi)
	559	IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
	560
	561	CALL gather_omp(phy_fsic, phy_mpi)
	562	IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
	563
	564	CALL gather_omp(phy_fter, phy_mpi)
	565	IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
	566
	567	CALL gather_omp(phy_flic, phy_mpi)
	568	IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
	569
	570	CALL gather_omp(phy_sst, phy_mpi)
	571	IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
	572
	573	CALL gather_omp(phy_bil, phy_mpi)
	574	IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
	575
	576	CALL gather_omp(phy_alb, phy_mpi)
	577	IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
	578
	579	CALL gather_omp(phy_rug, phy_mpi)
	580	IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
[4619]	581
[3927]	582	END SUBROUTINE writelim_unstruct
	583
	584
	585
	586	SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
	587	phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
	588
[5110]	589	USE lmdz_phys_para, ONLY: is_master
	590	USE lmdz_grid_phy, ONLY: klon_glo
	591	USE lmdz_phys_transfert_para, ONLY: gather
[3927]	592	USE phys_cal_mod, ONLY: year_len
[5088]	593	USE netcdf, ONLY:nf90_clobber,nf90_close,nf90_noerr,nf90_strerror,nf90_put_att,nf90_def_var,&
	594	nf90_def_dim,nf90_create,nf90_put_var,nf90_unlimited,nf90_global,nf90_64bit_offset,&
	595	nf90_enddef
[5090]	596	USE lmdz_cppkeys_wrapper, ONLY: nf90_format
[3927]	597	IMPLICIT NONE
	598
	599	INTEGER, INTENT (IN) :: klon
	600	REAL, INTENT (IN) :: phy_nat(klon, year_len)
	601	REAL, INTENT (IN) :: phy_alb(klon, year_len)
	602	REAL, INTENT (IN) :: phy_sst(klon, year_len)
	603	REAL, INTENT (IN) :: phy_bil(klon, year_len)
	604	REAL, INTENT (IN) :: phy_rug(klon, year_len)
	605	REAL, INTENT (IN) :: phy_ice(klon, year_len)
	606	REAL, INTENT (IN) :: phy_fter(klon, year_len)
	607	REAL, INTENT (IN) :: phy_foce(klon, year_len)
	608	REAL, INTENT (IN) :: phy_flic(klon, year_len)
	609	REAL, INTENT (IN) :: phy_fsic(klon, year_len)
	610
	611	REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
	612	! on the whole physics grid
	613	INTEGER :: k
	614	INTEGER ierr
	615	INTEGER dimfirst(3)
	616	INTEGER dimlast(3)
	617
	618	INTEGER nid, ndim, ntim
	619	INTEGER dims(2), debut(2), epais(2)
	620	INTEGER id_tim
	621	INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
	622	INTEGER id_fter, id_foce, id_fsic, id_flic
	623
	624	IF (is_master) THEN
	625
	626	PRINT *, 'writelim: Ecriture du fichier limit'
	627
[5088]	628	ierr = nf90_create('limit.nc', IOR(nf90_clobber,nf90_64bit_offset), nid)
[3927]	629
[5088]	630	ierr = nf90_put_att(nid, nf90_global, 'title', &
[3927]	631	'Fichier conditions aux limites')
[5113]	632	! ierr = nf90_def_dim (nid, "points_physiques", klon, ndim)
[5088]	633	ierr = nf90_def_dim(nid, 'points_physiques', klon_glo, ndim)
	634	ierr = nf90_def_dim(nid, 'time', nf90_unlimited, ntim)
[3927]	635
	636	dims(1) = ndim
	637	dims(2) = ntim
	638
[5088]	639	ierr = nf90_def_var(nid, 'TEMPS', nf90_format, ntim, id_tim)
	640	ierr = nf90_put_att(nid, id_tim, 'title', 'Jour dans l annee')
[3927]	641
[5088]	642	ierr = nf90_def_var(nid, 'NAT', nf90_format, dims, id_nat)
	643	ierr = nf90_put_att(nid, id_nat, 'title', &
[3927]	644	'Nature du sol (0,1,2,3)')
	645
[5088]	646	ierr = nf90_def_var(nid, 'SST', nf90_format, dims, id_sst)
	647	ierr = nf90_put_att(nid, id_sst, 'title', &
[3927]	648	'Temperature superficielle de la mer')
	649
[5088]	650	ierr = nf90_def_var(nid, 'BILS', nf90_format, dims, id_bils)
	651	ierr = nf90_put_att(nid, id_bils, 'title', &
[3927]	652	'Reference flux de chaleur au sol')
	653
[5088]	654	ierr = nf90_def_var(nid, 'ALB', nf90_format, dims, id_alb)
	655	ierr = nf90_put_att(nid, id_alb, 'title', 'Albedo a la surface')
[3927]	656
[5088]	657	ierr = nf90_def_var(nid, 'RUG', nf90_format, dims, id_rug)
	658	ierr = nf90_put_att(nid, id_rug, 'title', 'Rugosite')
[3927]	659
[5088]	660	ierr = nf90_def_var(nid, 'FTER', nf90_format, dims, id_fter)
	661	ierr = nf90_put_att(nid, id_fter, 'title','Frac. Land')
	662	ierr = nf90_def_var(nid, 'FOCE', nf90_format, dims, id_foce)
	663	ierr = nf90_put_att(nid, id_foce, 'title','Frac. Ocean')
	664	ierr = nf90_def_var(nid, 'FSIC', nf90_format, dims, id_fsic)
	665	ierr = nf90_put_att(nid, id_fsic, 'title','Frac. Sea Ice')
	666	ierr = nf90_def_var(nid, 'FLIC', nf90_format, dims, id_flic)
	667	ierr = nf90_put_att(nid, id_flic, 'title','Frac. Land Ice')
[3927]	668
[5088]	669	ierr = nf90_enddef(nid)
	670	IF (ierr/=nf90_noerr) THEN
[3927]	671	WRITE (, ) 'writelim error: failed to end define mode'
[5088]	672	WRITE (, ) nf90_strerror(ierr)
[3927]	673	END IF
	674
	675
	676	! write the 'times'
	677	DO k = 1, year_len
[5076]	678	ierr = nf90_put_var(nid, id_tim, k, [k])
[5088]	679	IF (ierr/=nf90_noerr) THEN
[3927]	680	WRITE (, ) 'writelim error with temps(k),k=', k
[5088]	681	WRITE (, ) nf90_strerror(ierr)
[3927]	682	END IF
	683	END DO
	684
	685	END IF ! of if (is_master)
	686
	687	! write the fields, after having collected them on master
	688
	689	CALL gather(phy_nat, phy_glo)
	690	IF (is_master) THEN
[5073]	691	ierr = nf90_put_var(nid, id_nat, phy_glo)
[5088]	692	IF (ierr/=nf90_noerr) THEN
[3927]	693	WRITE (, ) 'writelim error with phy_nat'
[5088]	694	WRITE (, ) nf90_strerror(ierr)
[3927]	695	END IF
	696	END IF
	697
	698	CALL gather(phy_sst, phy_glo)
	699	IF (is_master) THEN
[5073]	700	ierr = nf90_put_var(nid, id_sst, phy_glo)
[5088]	701	IF (ierr/=nf90_noerr) THEN
[3927]	702	WRITE (, ) 'writelim error with phy_sst'
[5088]	703	WRITE (, ) nf90_strerror(ierr)
[3927]	704	END IF
	705	END IF
	706
	707	CALL gather(phy_bil, phy_glo)
	708	IF (is_master) THEN
[5073]	709	ierr = nf90_put_var(nid, id_bils, phy_glo)
[5088]	710	IF (ierr/=nf90_noerr) THEN
[3927]	711	WRITE (, ) 'writelim error with phy_bil'
[5088]	712	WRITE (, ) nf90_strerror(ierr)
[3927]	713	END IF
	714	END IF
	715
	716	CALL gather(phy_alb, phy_glo)
	717	IF (is_master) THEN
[5073]	718	ierr = nf90_put_var(nid, id_alb, phy_glo)
[5088]	719	IF (ierr/=nf90_noerr) THEN
[3927]	720	WRITE (, ) 'writelim error with phy_alb'
[5088]	721	WRITE (, ) nf90_strerror(ierr)
[3927]	722	END IF
	723	END IF
	724
	725	CALL gather(phy_rug, phy_glo)
	726	IF (is_master) THEN
[5073]	727	ierr = nf90_put_var(nid, id_rug, phy_glo)
[5088]	728	IF (ierr/=nf90_noerr) THEN
[3927]	729	WRITE (, ) 'writelim error with phy_rug'
[5088]	730	WRITE (, ) nf90_strerror(ierr)
[3927]	731	END IF
	732	END IF
	733
	734	CALL gather(phy_fter, phy_glo)
	735	IF (is_master) THEN
[5073]	736	ierr = nf90_put_var(nid, id_fter, phy_glo)
[5088]	737	IF (ierr/=nf90_noerr) THEN
[3927]	738	WRITE (, ) 'writelim error with phy_fter'
[5088]	739	WRITE (, ) nf90_strerror(ierr)
[3927]	740	END IF
	741	END IF
	742
	743	CALL gather(phy_foce, phy_glo)
	744	IF (is_master) THEN
[5073]	745	ierr = nf90_put_var(nid, id_foce, phy_glo)
[5088]	746	IF (ierr/=nf90_noerr) THEN
[3927]	747	WRITE (, ) 'writelim error with phy_foce'
[5088]	748	WRITE (, ) nf90_strerror(ierr)
[3927]	749	END IF
	750	END IF
	751
	752	CALL gather(phy_fsic, phy_glo)
	753	IF (is_master) THEN
[5073]	754	ierr = nf90_put_var(nid, id_fsic, phy_glo)
[5088]	755	IF (ierr/=nf90_noerr) THEN
[3927]	756	WRITE (, ) 'writelim error with phy_fsic'
[5088]	757	WRITE (, ) nf90_strerror(ierr)
[3927]	758	END IF
	759	END IF
	760
	761	CALL gather(phy_flic, phy_glo)
	762	IF (is_master) THEN
[5073]	763	ierr = nf90_put_var(nid, id_flic, phy_glo)
[5088]	764	IF (ierr/=nf90_noerr) THEN
[3927]	765	WRITE (, ) 'writelim error with phy_flic'
[5088]	766	WRITE (, ) nf90_strerror(ierr)
[3927]	767	END IF
	768	END IF
	769
	770	! close file:
	771	IF (is_master) THEN
[5088]	772	ierr = nf90_close(nid)
[3927]	773	END IF
	774
	775	END SUBROUTINE writelim
	776
	777	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	778
	779	SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
	780	USE dimphy
	781	USE phys_cal_mod , ONLY: year_len
	782	IMPLICIT NONE
	783
	784	INTEGER nlon, type_profil, i, k, j
	785	REAL :: rlatd(nlon), phy_sst(nlon, year_len)
	786	INTEGER imn, imx, amn, amx, kmn, kmx
	787	INTEGER p, pplus, nlat_max
	788	PARAMETER (nlat_max=72)
	789	REAL x_anom_sst(nlat_max)
	790	CHARACTER (LEN=20) :: modname='profil_sst'
	791	CHARACTER (LEN=80) :: abort_message
	792
	793	IF (klon/=nlon) THEN
	794	abort_message='probleme de dimensions dans profil_sst'
	795	CALL abort_physic(modname,abort_message,1)
	796	ENDIF
	797	WRITE (, ) ' profil_sst: type_profil=', type_profil
	798	DO i = 1, year_len
	799	! phy_sst(:,i) = 260.+50.cos(rlatd(:))*2
	800
	801	! Rajout fbrlmd
	802
	803	IF (type_profil==1) THEN
	804	! Méthode 1 "Control" faible plateau à l'Equateur
	805	DO j = 1, klon
	806	phy_sst(j, i) = 273. + 27.(1-sin(1.5rlatd(j))**2)
	807	! PI/3=1.047197551
	808	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	809	phy_sst(j, i) = 273.
	810	END IF
	811	END DO
	812	END IF
	813	IF (type_profil==2) THEN
	814	! Méthode 2 "Flat" fort plateau à l'Equateur
	815	DO j = 1, klon
	816	phy_sst(j, i) = 273. + 27.(1-sin(1.5rlatd(j))**4)
	817	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	818	phy_sst(j, i) = 273.
	819	END IF
	820	END DO
	821	END IF
	822
	823
	824	IF (type_profil==3) THEN
	825	! Méthode 3 "Qobs" plateau réel à l'Equateur
	826	DO j = 1, klon
	827	phy_sst(j, i) = 273. + 0.527.(2-sin(1.5rlatd(j))2-sin(1.5 &
	828	rlatd(j))**4)
	829	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	830	phy_sst(j, i) = 273.
	831	END IF
	832	END DO
	833	END IF
	834
	835	IF (type_profil==4) THEN
	836	! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
	837	DO j = 1, klon
	838	phy_sst(j, i) = 273. + 0.529.(2-sin(1.5rlatd(j))2-sin(1.5 &
	839	rlatd(j))**4)
	840	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	841	phy_sst(j, i) = 273.
	842	END IF
	843	END DO
	844	END IF
	845
	846	IF (type_profil==5) THEN
	847	! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
	848	DO j = 1, klon
	849	phy_sst(j, i) = 273. + 2. + 0.527.(2-sin(1.5rlatd(j))*2-sin(1.5 &
	850	rlatd(j))*4)
	851	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	852	phy_sst(j, i) = 273. + 2.
	853	END IF
	854
	855	END DO
	856	END IF
	857
	858	IF (type_profil==6) THEN
	859	! Méthode 6 "cst" valeur constante de SST
	860	DO j = 1, klon
	861	phy_sst(j, i) = 288.
	862	END DO
	863	END IF
	864
	865
	866	IF (type_profil==7) THEN
	867	! Méthode 7 "cst" valeur constante de SST +2
	868	DO j = 1, klon
	869	phy_sst(j, i) = 288. + 2.
	870	END DO
	871	END IF
	872
	873	p = 0
	874	IF (type_profil==8) THEN
	875	! Méthode 8 profil anomalies SST du modèle couplé AR4
	876	DO j = 1, klon
	877	IF (rlatd(j)==rlatd(j-1)) THEN
	878	phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
	879	0.527.(2-sin(1.5rlatd(j))2-sin(1.5rlatd(j))**4)
	880	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	881	phy_sst(j, i) = 273. + x_anom_sst(pplus)
	882	END IF
	883	ELSE
	884	p = p + 1
	885	pplus = 73 - p
	886	phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
	887	0.527.(2-sin(1.5rlatd(j))2-sin(1.5rlatd(j))**4)
	888	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	889	phy_sst(j, i) = 273. + x_anom_sst(pplus)
	890	END IF
	891	WRITE (, ) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
	892	END IF
	893	END DO
	894	END IF
	895
	896	IF (type_profil==9) THEN
	897	! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
	898	DO j = 1, klon
	899	phy_sst(j, i) = 273. - 2. + 0.527.(2-sin(1.5rlatd(j))*2-sin(1.5 &
	900	rlatd(j))*4)
	901	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	902	phy_sst(j, i) = 273. - 2.
	903	END IF
	904	END DO
	905	END IF
	906
	907
	908	IF (type_profil==10) THEN
	909	! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
	910	DO j = 1, klon
	911	phy_sst(j, i) = 273. + 4. + 0.527.(2-sin(1.5rlatd(j))*2-sin(1.5 &
	912	rlatd(j))*4)
	913	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	914	phy_sst(j, i) = 273. + 4.
	915	END IF
	916	END DO
	917	END IF
	918
	919	IF (type_profil==11) THEN
	920	! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
	921	DO j = 1, klon
	922	phy_sst(j, i) = 273. + 0.527.(2-sin(1.5rlatd(j))2-sin(1.5 &
	923	rlatd(j))**4)
	924	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	925	phy_sst(j, i) = 273.
	926	END IF
	927	END DO
	928	END IF
	929
	930	IF (type_profil==12) THEN
	931	! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
	932	DO j = 1, klon
	933	phy_sst(j, i) = 273. + 4. + 0.527.(2-sin(1.5rlatd(j))*2-sin(1.5 &
	934	rlatd(j))*4)
	935	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	936	phy_sst(j, i) = 273. + 4.
	937	END IF
	938	END DO
	939	END IF
	940
	941	IF (type_profil==13) THEN
	942	! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
	943	DO j = 1, klon
	944	phy_sst(j, i) = 273. + 0.529.(2-sin(1.5rlatd(j))2-sin(1.5 &
	945	rlatd(j))**4)
	946	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	947	phy_sst(j, i) = 273.
	948	END IF
	949	END DO
	950	END IF
	951
	952	IF (type_profil==14) THEN
	953	! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
	954	DO j = 1, klon
	955	phy_sst(j, i) = 273. + 2. + 0.529.(2-sin(1.5rlatd(j))*2-sin(1.5 &
	956	rlatd(j))*4)
	957	IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
	958	phy_sst(j, i) = 273.
	959	END IF
	960	END DO
	961	END IF
	962
[5117]	963	IF (type_profil==20) THEN
[5103]	964	PRINT*,'Profile SST 20'
[5093]	965	! Méthode 13 "Qmax2K" plateau réel é\| l'Equateur augmenté +2K
[3927]	966
	967	do j=1,klon
	968	phy_sst(j,i)=248.+55.(1-sin(rlatd(j))*2)
	969	enddo
	970	endif
	971
[5117]	972	IF (type_profil==21) THEN
[5103]	973	PRINT*,'Profile SST 21'
[5093]	974	! Méthode 13 "Qmax2K" plateau réel é\| l'Equateur augmenté +2K
[3927]	975	do j=1,klon
	976	phy_sst(j,i)=252.+55.(1-sin(rlatd(j))*2)
	977	enddo
	978	endif
	979
	980
	981
	982	END DO
	983
	984	! IM beg : verif profil SST: phy_sst
	985	amn = min(phy_sst(1,1), 1000.)
	986	amx = max(phy_sst(1,1), -1000.)
	987	imn = 1
	988	kmn = 1
	989	imx = 1
	990	kmx = 1
	991	DO k = 1, year_len
	992	DO i = 2, nlon
	993	IF (phy_sst(i,k)<amn) THEN
	994	amn = phy_sst(i, k)
	995	imn = i
	996	kmn = k
	997	END IF
	998	IF (phy_sst(i,k)>amx) THEN
	999	amx = phy_sst(i, k)
	1000	imx = i
	1001	kmx = k
	1002	END IF
	1003	END DO
	1004	END DO
	1005
	1006	PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
	1007	PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
	1008	! IM end : verif profil SST: phy_sst
	1009
[5105]	1010
[3927]	1011	END SUBROUTINE profil_sst
	1012
	1013	END MODULE phyaqua_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format