Context Navigation

phyaqua.F @ 1529

Last change on this file since 1529 was 1529, checked in by Laurent Fairhead, 13 years ago

Necessary modifications to use the model in an aquaplanet or terraplanet configuration

Modifications nécessaires pour l'utilisation du moèle en configuration aquaplanète
ou terraplanète

Dans la dynamique:
Changement pour le paramètre iflag_phys:

auparavant : 0 pas de physique, 1 phyique, 2 rappel
ici on ajoute iflag_phys>=100 pour les aqua et terra

Du coup, dans leapfrog.F et gcm.F on appelle la physique si
iflag_phys=1.or.iflag_phys>=100
Dans iniacademic, on initialise si iflag_phys>=2 au lieu de =2
Dans gcm.F, on appelle en plus iniaqua (sous une clef CPP_EARTH)
Dans iniacademic, on met ps=108080 pour les aqua et terra pour répondre
à une specification CFMIP.

Dans la physique:
On ajoute phyaqua.F qui contient :
iniaqua : initialise les startphy.nc et limit.nc pour la phyique
zenang_an : calcule un ensolleillement moyen sur l'année en fonction de

la latiude (A. Campoy).

profil_sst : qui calcule différents profils latitudinaux de SSTs.
writelim : écrit le fichier limit.nc

Dans physiq.F
On ajoute la possibilité d'appeller un calcul d'ensoleillement moyen
sur l'année quand solarlong0=1000.

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

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

Download in other formats:

Original Format