source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/phyaqua.F @ 2525

Last change on this file since 2525 was 1785, checked in by Ehouarn Millour, 11 years ago

Transformation de l'include indicesol.h en un module indice_sol_mod et modification des appels dans tous les fichiers concernés.
Aucun changement des résultats ni des sorties du modèle vs 1784.
UG

...................................................

Replacement of the indicesol.h include by a module named indice_sol_mod. Modification of the calls in every affected files.
Results and outputs of simulations are unchanged in comparison with rev 1784.
UG
File size: 27.1 KB
Line 
1! Routines complementaires pour la physique planetaire.
2
3
4      subroutine iniaqua(nlon,latfi,lonfi,iflag_phys)
5
6!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
7!  Creation d'un etat initial et de conditions aux limites
8!  (resp startphy.nc et limit.nc) pour des configurations idealisees
9! du modele LMDZ dans sa version terrestre.
10!  iflag_phys est un parametre qui controle
11!  iflag_phys = N 
12!    de 100 a 199 : aqua planetes avec SST forcees
13!                 N-100 determine le type de SSTs
14!    de 200 a 299 : terra planetes avec Ts calcule
15!       
16!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
17
18      USE comgeomphy, 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      USE indice_sol_mod
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 "temps.h"
37
38      integer,intent(in) :: nlon,iflag_phys
39cIM ajout latfi, lonfi
40      real,intent(in) :: lonfi(nlon),latfi(nlon)
41
42      INTEGER type_profil,type_aqua
43
44c  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
55c    Arguments :
56c    -----------
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
114c-----------------------------------------------------------------------
115c   dynamial tendencies :
116c   ---------------------
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
151cIM 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!---------------------------------------------------------------------
174c Creation des conditions aux limites:
175c ------------------------------------
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
195c  Conditions aux limites:
196c  -----------------------
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
245cIM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
246cIM 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
257cIM 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!---------------------------------------------------------------------
266c Ecriture de l'etat initial:
267c ---------------------------
268
269C
270C Ecriture etat initial physique
271C
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 )
280c     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
289c=======================================================================
290c  Profils initiaux
291c=======================================================================
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
364c====================================================================
365c====================================================================
366      SUBROUTINE zenang_an(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
367      USE dimphy
368      IMPLICIT none
369c====================================================================
370c=============================================================
371c         CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
372c Auteur : A. Campoy et F. Hourdin
373c Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
374c          et l'ensoleillement moyen entre gmtime1 et gmtime2
375c          connaissant la declinaison, la latitude et la longitude.
376c
377c   Dans cette version particuliere, on calcule le rayonnement
378c  moyen sur l'année à chaque latitude.
379c   angle zenithal calculé pour obtenir un
380c   Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
381c   en moyenne annuelle.
382c   Spécifique de la terre. Utilisé pour les aqua planetes.
383c
384c Rque   : Different de la routine angle en ce sens que zenang
385c          fournit des moyennes de pmu0 et non des valeurs
386c          instantanees, du coup frac prend toutes les valeurs
387c          entre 0 et 1.
388c Date   : premiere version le 13 decembre 1994
389c          revu pour  GCM  le 30 septembre 1996
390c===============================================================
391c longi----INPUT : la longitude vraie de la terre dans son plan
392c                  solaire a partir de l'equinoxe de printemps (degre)
393c gmtime---INPUT : temps universel en fraction de jour
394c pdtrad---INPUT : pas de temps du rayonnement (secondes)
395c lat------INPUT : latitude en degres
396c long-----INPUT : longitude en degres
397c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
398c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
399c================================================================
400#include "YOMCST.h"
401c================================================================
402      logical cycle_diurne
403      real  gmtime
404      real rlat(klon), rlon(klon), rmu0(klon), fract(klon)
405c================================================================
406      integer i
407      real gmtime1, gmtime2
408      real pi_local
409
410
411      real rmu0m(klon),rmu0a(klon)
412c
413
414      pi_local = 4.0 * ATAN(1.0)
415
416c================================================================
417c  Calcul de l'angle zenithal moyen sur la journee
418c================================================================
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
434c================================================================
435!  Avec ou sans cycle diurne
436c================================================================
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)
484c
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)
509c
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'
519c
520        ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
521c
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)
527c
528        dims(1) = ndim
529        dims(2) = ntim
530c
531ccc      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")
535ccc      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)")
539ccc      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")
543ccc      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")
547ccc      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")
551ccc      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")
564c
565        ierr = NF_ENDDEF(nid)
566c
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
710c
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)
717c
718      if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua'
719      do i=1,360
720c      phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
721
722c Rajout fbrlmd
723
724      if(type_profil.EQ.1)then
725c     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)
728c        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
735c     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
746c     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
757c     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
768c     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
780c     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
788c     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
796c     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
818c     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
830c     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
841c     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
852c     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
863c     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
874c     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
886cIM 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
903c
904       PRINT*,' debut avant writelim min max phy_sst',imn,kmn,amn,
905     & imx,kmx,amx
906cIM end : verif profil SST: phy_sst
907
908       return
909       end
Note: See TracBrowser for help on using the repository browser.