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