source: LMDZ5/branches/LF-private/libf/phylmd/phyaqua_mod.F @ 4400

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

OpenMP bug fix in iniaqua: local variables are shared by all threads only if they are static (ie with "save" attribute).
While at it, changed phyaqua.F into module phyaqua_mod.F
EM

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