source: LMDZ5/trunk/libf/phylmd/phyaqua.F @ 1544

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

Some initialisations missing for the aquaplanet configuration

Initialisations manquantes pour les aquaplanètes
File size: 23.6 KB
RevLine 
[1529]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
39cIM ajout latfi, lonfi
40      REAL, DIMENSION (nlon) :: lonfi, latfi
41      INTEGER type_profil,type_aqua
42
43c  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
54c    Arguments :
55c    -----------
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
111c-----------------------------------------------------------------------
112c   dynamial tendencies :
113c   ---------------------
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
145cIM 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!---------------------------------------------------------------------
163c Creation des conditions aux limites:
164c ------------------------------------
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
177c  Conditions aux limites:
178c  -----------------------
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
219cIM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
220cIM 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
231cIM 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!---------------------------------------------------------------------
240c Ecriture de l'etat initial:
241c ---------------------------
242
243C
244C Ecriture etat initial physique
245C
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 )
254c     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
263c=======================================================================
264c  Profils initiaux
265c=======================================================================
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
[1530]312      falb1=albedo
313      falb2=albedo
314      zmax0=0.
315      f0=0.
316      ema_work1=0.
317      ema_work2=0.
318      wake_deltat=0.
319      wake_deltaq=0.
320      wake_s=0.
321      wake_cstar=0.
322      wake_pe=0.
323      wake_fip=0.
324      fm_therm=0.
325      entr_therm=0.
326      detr_therm=0.
327
328
[1529]329      CALL phyredem ("startphy.nc")
330
331        print*,'apres phyredem'
332      call phys_state_var_end
333
334      return
335      end
336
337
338c====================================================================
339c====================================================================
340      SUBROUTINE zenang_an(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
341      USE dimphy
342      IMPLICIT none
343c====================================================================
344c=============================================================
345c         CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
346c Auteur : A. Campoy et F. Hourdin
347c Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
348c          et l'ensoleillement moyen entre gmtime1 et gmtime2
349c          connaissant la declinaison, la latitude et la longitude.
350c
351c   Dans cette version particuliere, on calcule le rayonnement
352c  moyen sur l'année à chaque latitude.
353c   angle zenithal calculé pour obtenir un
354c   Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
355c   en moyenne annuelle.
356c   Spécifique de la terre. Utilisé pour les aqua planetes.
357c
358c Rque   : Different de la routine angle en ce sens que zenang
359c          fournit des moyennes de pmu0 et non des valeurs
360c          instantanees, du coup frac prend toutes les valeurs
361c          entre 0 et 1.
362c Date   : premiere version le 13 decembre 1994
363c          revu pour  GCM  le 30 septembre 1996
364c===============================================================
365c longi----INPUT : la longitude vraie de la terre dans son plan
366c                  solaire a partir de l'equinoxe de printemps (degre)
367c gmtime---INPUT : temps universel en fraction de jour
368c pdtrad---INPUT : pas de temps du rayonnement (secondes)
369c lat------INPUT : latitude en degres
370c long-----INPUT : longitude en degres
371c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
372c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
373c================================================================
374#include "YOMCST.h"
375c================================================================
376      logical cycle_diurne
377      real  gmtime
378      real rlat(klon), rlon(klon), rmu0(klon), fract(klon)
379c================================================================
380      integer i
381      real gmtime1, gmtime2
382      real pi_local
383
384
385      real rmu0m(klon),rmu0a(klon)
386c
387
388      pi_local = 4.0 * ATAN(1.0)
389
390c================================================================
391c  Calcul de l'angle zenithal moyen sur la journee
392c================================================================
393
394      DO i=1,klon
395        fract(i)=1.
396!  Calcule du flux moyen
397        IF (abs(rlat(i)).LE.28.75) THEN
398           rmu0m(i)=(210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
399        ELSEIF (abs(rlat(i)).LE.43.75) THEN
400          rmu0m(i)=(187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
401        ELSEIF (abs(rlat(i)).LE.71.25) THEN
402          rmu0m(i)=(162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
403        ELSE
404          rmu0m(i)=(172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
405        ENDIF
406      ENDDO
407
408c================================================================
409!  Avec ou sans cycle diurne
410c================================================================
411
412      IF (cycle_diurne) THEN
413
414!  On redecompose flux  au sommet suivant un cycle diurne idealise
415!  identique a toutes les latitudes.
416
417         DO i=1,klon
418           rmu0a(i)=2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
419           rmu0(i)=rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*
420     &      rlon(i)/360.))-rmu0a(i)/sqrt(2.)
421         ENDDO
422
423         DO i=1,klon
424           IF (rmu0(i).LE.0.) THEN
425              rmu0(i)=0.
426              fract(i)=0.
427           ELSE
428              fract(i)=1.
429           ENDIF
430         ENDDO
431
432!  Affichage de l'angel zenitale
433!     print*,'************************************'
434!     print*,'************************************'
435!     print*,'************************************'
436!     print*,'latitude=',rlat(i),'longitude=',rlon(i)
437!     print*,'rmu0m=',rmu0m(i)
438!     print*,'rmu0a=',rmu0a(i)
439!     print*,'rmu0=',rmu0(i)
440                                                               
441      ELSE
442
443        DO i=1,klon
444           fract(i)=0.5
445           rmu0(i)=rmu0m(i)*2.
446        ENDDO
447
448      ENDIF
449
450      RETURN
451      END
452      subroutine writelim
453     s   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
454     s    phy_fter,phy_foce,phy_flic,phy_fsic)
455c
456!#include "dimensions.h"
457!#include "dimphy.h"
458#include "netcdf.inc"
459 
460      integer klon
461      REAL phy_nat(klon,360)
462      REAL phy_alb(klon,360)
463      REAL phy_sst(klon,360)
464      REAL phy_bil(klon,360)
465      REAL phy_rug(klon,360)
466      REAL phy_ice(klon,360)
467      REAL phy_fter(klon,360)
468      REAL phy_foce(klon,360)
469      REAL phy_flic(klon,360)
470      REAL phy_fsic(klon,360)
471 
472      INTEGER ierr
473      INTEGER dimfirst(3)
474      INTEGER dimlast(3)
475c
476      INTEGER nid, ndim, ntim
477      INTEGER dims(2), debut(2), epais(2)
478      INTEGER id_tim
479      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
480      INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
481 
482      PRINT*, 'Ecriture du fichier limit'
483c
484      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
485c
486      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
487     .                       "Fichier conditions aux limites")
488      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
489      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
490c
491      dims(1) = ndim
492      dims(2) = ntim
493c
494ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
495      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
496      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
497     .                        "Jour dans l annee")
498ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
499      ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
500      ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
501     .                        "Nature du sol (0,1,2,3)")
502ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
503      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
504      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
505     .                        "Temperature superficielle de la mer")
506ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
507      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
508      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
509     .                        "Reference flux de chaleur au sol")
510ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
511      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
512      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
513     .                        "Albedo a la surface")
514ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
515      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
516      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
517     .                        "Rugosite")
518
519      ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
520      ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
521      ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
522      ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
523      ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
524      ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
525      ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
526      ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
527c
528      ierr = NF_ENDDEF(nid)
529c
530      DO k = 1, 360
531c
532      debut(1) = 1
533      debut(2) = k
534      epais(1) = klon
535      epais(2) = 1
536c
537      print*,'Instant ',k
538#ifdef NC_DOUBLE
539      print*,'NC DOUBLE'
540      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
541      ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
542      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
543      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
544      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
545      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
546      ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais,phy_fter(1,k))
547      ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais,phy_foce(1,k))
548      ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais,phy_fsic(1,k))
549      ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais,phy_flic(1,k))
550#else
551      print*,'NC PAS DOUBLE'
552      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
553      ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
554      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
555      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
556      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
557      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
558      ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais,phy_fter(1,k))
559      ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais,phy_foce(1,k))
560      ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais,phy_fsic(1,k))
561      ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais,phy_flic(1,k))
562
563#endif
564c
565      ENDDO
566c
567      ierr = NF_CLOSE(nid)
568c
569      return
570      end
571
572      SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
573      use dimphy
574      IMPLICIT none
575c
576      INTEGER nlon, type_profil, i, k, j
577      REAL :: rlatd(nlon), phy_sst(nlon, 360)
578      INTEGER imn, imx, amn, amx, kmn, kmx
579      INTEGER p, pplus, nlat_max
580      parameter (nlat_max=72)
581      REAL x_anom_sst(nlat_max)
582c
583      if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua'
584      do i=1,360
585c      phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
586
587c Rajout fbrlmd
588
589      if(type_profil.EQ.1)then
590c     Méthode 1 "Control" faible plateau à l'Equateur
591      do j=1,klon
592       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**2)
593c        PI/3=1.047197551
594      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
595         phy_sst(j,i)=273.
596        endif
597      enddo
598      endif
599      if(type_profil.EQ.2)then
600c     Méthode 2 "Flat" fort plateau à l'Equateur
601      do j=1,klon
602       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**4)
603      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
604         phy_sst(j,i)=273.
605        endif
606      enddo
607      endif
608
609
610      if (type_profil.EQ.3) then
611c     Méthode 3 "Qobs" plateau réel à l'Equateur
612      do j=1,klon
613        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2
614     &                 -sin(1.5*rlatd(j))**4)
615      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
616         phy_sst(j,i)=273.
617        endif
618      enddo
619      endif
620
621      if (type_profil.EQ.4) then
622c     Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
623      do j=1,klon
624        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2
625     &                 -sin(1.5*rlatd(j))**4)
626      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
627         phy_sst(j,i)=273.
628        endif
629      enddo
630      endif
631
632      if (type_profil.EQ.5) then
633c     Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
634      do j=1,klon
635        phy_sst(j,i)=273.+2.+0.5*27.*(2-sin(1.5*rlatd(j))**2
636     &                 -sin(1.5*rlatd(j))**4)
637      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
638         phy_sst(j,i)=273.+2.
639        endif
640
641      enddo
642      endif
643
644      if(type_profil.EQ.6)then
645c     Méthode 6 "cst" valeur constante de SST
646      do j=1,klon
647       phy_sst(j,i)=288.
648      enddo
649      endif
650
651
652        if(type_profil.EQ.7)then
653c     Méthode 7 "cst" valeur constante de SST +2
654      do j=1,klon
655       phy_sst(j,i)=288.+2.
656      enddo
657      endif
658
659        p=0
660        if(type_profil.EQ.8)then
661c     Méthode 8 profil anomalies SST du modèle couplé AR4
662       do j=1,klon
663         if (rlatd(j).EQ.rlatd(j-1)) then
664       phy_sst(j,i)=273.+x_anom_sst(pplus)
665     &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
666          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
667            phy_sst(j,i)=273.+x_anom_sst(pplus)
668          endif
669        else
670          p=p+1
671          pplus=73-p
672        phy_sst(j,i)=273.+x_anom_sst(pplus)
673     &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
674          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
675            phy_sst(j,i)=273.+x_anom_sst(pplus)
676          endif
677          write (*,*) rlatd(j),x_anom_sst(pplus),phy_sst(j,i)
678        endif
679      enddo
680      endif
681
682      if (type_profil.EQ.9) then
683c     Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
684      do j=1,klon
685        phy_sst(j,i)=273.-2.+0.5*27.*(2-sin(1.5*rlatd(j))**2
686     &                 -sin(1.5*rlatd(j))**4)
687      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
688         phy_sst(j,i)=273.-2.
689        endif
690      enddo
691      endif
692
693
694      if (type_profil.EQ.10) then
695c     Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
696      do j=1,klon
697        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2
698     &                 -sin(1.5*rlatd(j))**4)
699        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
700         phy_sst(j,i)=273.+4.
701        endif
702      enddo
703      endif
704
705      if (type_profil.EQ.11) then
706c     Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
707      do j=1,klon
708        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2
709     &                 -sin(1.5*rlatd(j))**4)
710        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
711         phy_sst(j,i)=273.
712        endif
713      enddo
714      endif
715
716      if (type_profil.EQ.12) then
717c     Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
718      do j=1,klon
719        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2
720     &                 -sin(1.5*rlatd(j))**4)
721        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
722         phy_sst(j,i)=273.+4.
723        endif
724      enddo
725      endif
726
727      if (type_profil.EQ.13) then
728c     Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
729      do j=1,klon
730        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2
731     &                 -sin(1.5*rlatd(j))**4)
732      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
733         phy_sst(j,i)=273.
734        endif
735      enddo
736      endif
737
738      if (type_profil.EQ.14) then
739c     Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
740      do j=1,klon
741        phy_sst(j,i)=273.+2.+0.5*29.*(2-sin(1.5*rlatd(j))**2
742     &                 -sin(1.5*rlatd(j))**4)
743      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
744         phy_sst(j,i)=273.
745        endif
746      enddo
747      endif
748
749      enddo
750
751cIM beg : verif profil SST: phy_sst
752       amn=MIN(phy_sst(1,1),1000.)
753       amx=MAX(phy_sst(1,1),-1000.)
754       DO k=1, 360
755       DO i=2, nlon
756        IF(phy_sst(i,k).LT.amn) THEN
757         amn=phy_sst(i,k)
758         imn=i
759         kmn=k
760        ENDIF
761        IF(phy_sst(i,k).GT.amx) THEN
762         amx=phy_sst(i,k)
763         imx=i
764         kmx=k
765        ENDIF
766       ENDDO
767       ENDDO
768c
769       PRINT*,' debut avant writelim min max phy_sst',imn,kmn,amn,
770     & imx,kmx,amx
771cIM end : verif profil SST: phy_sst
772
773       return
774       end
Note: See TracBrowser for help on using the repository browser.