Context Navigation

phyaqua_mod.F @ 1891

Last change on this file since 1891 was 1851, checked in by Ehouarn Millour, 11 years ago
OpenMP bug fix in iniaqua: local variables are shared by all threads only if they are static (ie with "save" attribute). While at it, changed phyaqua.F into module phyaqua_mod.F EM
File size: 28.0 KB

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

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

Download in other formats:

Original Format