Context Navigation

phyaqua_mod.F @ 3659

Last change on this file since 3659 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

Line
1	module phyaqua_mod
2	! Routines complementaires pour la physique planetaire.
3	implicit none
4
5	contains
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
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
25	USE fonte_neige_mod, only : fonte_neige_init
26	USE phys_state_var_mod
27	USE control_mod, only : dayref,nday,iphysiq
28	USE indice_sol_mod
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"
39	#include "temps.h"
40
41	integer,intent(in) :: nlon,iflag_phys
42	cIM ajout latfi, lonfi
43	real,intent(in) :: lonfi(nlon),latfi(nlon)
44
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
92	integer, save:: read_climoz=0 ! read ozone climatology
93
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
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
139	print*,'iniaqua:type_aqua, type_profil',type_aqua, type_profil
140
141	if (klon.ne.nlon) then
142	write(,)"iniaqua: klon=",klon," nlon=",nlon
143	stop'probleme de dimensions dans iniaqua'
144	endif
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
165	day_end=day_ini+nday
166	airefi=1.
167	zcufi=1.
168	zcvfi=1.
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
175
176	!---------------------------------------------------------------------
177	c Creation des conditions aux limites:
178	c ------------------------------------
179	! Initialisations des constantes
180	! Ajouter les manquants dans planete.def... (albedo etc)
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
192	write(,)"iniaqua: co2_ppm=",co2_ppm
193	solaire=solaire_omp
194	write(,)"iniaqua: solaire=",solaire
195	alb_ocean=alb_ocean_omp
196	write(,)"iniaqua: alb_ocean=",alb_ocean
197
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)
206	rugoro = 0.0
207	u_ancien = 0.0
208	v_ancien = 0.0
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
231	!$OMP MASTER
232	rugos_omp=rugos
233	CALL getin('rugos',rugos_omp)
234	!$OMP END MASTER
235	!$OMP BARRIER
236	rugos=rugos_omp
237	write(,) "iniaqua: rugos=",rugos
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
343	print*,'iniaqua: before phyredem'
344
345	falb1=albedo
346	falb2=albedo
347	zmax0=0.
348	f0=0.
349	sig1=0.
350	w01=0.
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
362	CALL phyredem ("startphy.nc")
363
364	print*,'iniaqua: after phyredem'
365	call phys_state_var_end
366
367	return
368	end subroutine iniaqua
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
484	END subroutine zenang_an
485
486	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
487
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
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
495	implicit none
496	!#include "dimensions.h"
497	!#include "dimphy.h"
498	#include "netcdf.inc"
499
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
514	integer :: k
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
525	if (is_mpi_root.and.is_omp_root) then
526
527	PRINT*, 'writelim: Ecriture du fichier limit'
528	c
529	ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
530	c
531	ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
532	. "Fichier conditions aux limites")
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)
536	c
537	dims(1) = ndim
538	dims(2) = ntim
539	c
540	ccc ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
541	ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
542	ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
543	. "Jour dans l annee")
544	ccc ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
545	ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
546	ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
547	. "Nature du sol (0,1,2,3)")
548	ccc ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
549	ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
550	ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
551	. "Temperature superficielle de la mer")
552	ccc ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
553	ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
554	ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
555	. "Reference flux de chaleur au sol")
556	ccc ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
557	ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
558	ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
559	. "Albedo a la surface")
560	ccc ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
561	ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
562	ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
563	. "Rugosite")
564
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")
573	c
574	ierr = NF_ENDDEF(nid)
575	if (ierr.ne.NF_NOERR) then
576	write(,) "writelim error: failed to end define mode"
577	write(,) NF_STRERROR(ierr)
578	endif
579	c
580
581	! write the 'times'
582	do k=1,360
583	#ifdef NC_DOUBLE
584	ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
585	#else
586	ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
587	#endif
588	if (ierr.ne.NF_NOERR) then
589	write(,) "writelim error with temps(k),k=",k
590	write(,) NF_STRERROR(ierr)
591	endif
592	enddo
593
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)
604	#endif
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
720	end subroutine writelim
721
722	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
723
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'
736	write(,)" profil_sst: type_profil=",type_profil
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.)
907	imn=1 ; kmn=1 ; imx=1 ; kmx=1
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
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
925	cIM end : verif profil SST: phy_sst
926
927	return
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