source: lmdz_wrf/branches/LMDZ_WRFmeas/WRFV3/lmdz/phyaqua.F90 @ 27

Last change on this file since 27 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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