Context Navigation

phyaqua.F @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur la r1706

Testing release based on r1706

File size: 26.6 KB

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

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

Download in other formats:

Original Format