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

Last change on this file since 1687 was 1671, checked in by Ehouarn Millour, 12 years ago
  • fixed "aquaplanet case" so that initializations (creation of files startphy.nc and limit.nc) now also works in parallel (mpi,omp,mixed).
  • call to "iniaqua" is now done from within "iniphysiq" ; also added some tests to check consistency between essential variables shared by dynamics and physics (planetary radius, gravity, Cp, ...)
  • simillarily adapted "phydev" routines, and added necessary routines to also be able to generate startphy.nc/restartphy.nc files there. Also removed common file "comcstphy.h" and replaced it with a module "comcstphy.F90"

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