Context Navigation

phyaqua_mod.F90 @ 5276

Last change on this file since 5276 was 5274, checked in by abarral, 23 hours ago
Replace yomcst.h by existing module
File size: 32.1 KB

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

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

Download in other formats:

Original Format